\documentclass[reprint,amsmath,amssymb,aps,prx,superscriptaddress]{revtex4-2}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{hyperref}

% \newcommand{\supp}{See Supplementary Material [url] for detailed proofs for mathematical conclusions in the main text, numerical results of quasi-1D and strip winding numbers, and the derivation of SGBZs for 2D HN model, which includes Ref.\cite{hatcherAlgebraicTopology2001}}

\begin{document}
\title{Non-Bloch Band Theory for 2D Geometry-Dependent Non-Hermitian Skin
Effect}
\author{Chenyang Wang}
\affiliation{State Key Laboratory of Low-Dimensional Quantum Physics, Department
of Physics, Tsinghua University, Beijing 100084, China}
\author{Jinghui Pi}
\affiliation{The Chinese University of Hong Kong Shenzhen Research Institute, 518057
Shenzhen, China}
\author{Qinxin Liu}
\affiliation{State Key Laboratory of Low-Dimensional Quantum Physics, Department
of Physics, Tsinghua University, Beijing 100084, China}
\author{Yaohua Li}
\affiliation{State Key Laboratory of Low-Dimensional Quantum Physics, Department
of Physics, Tsinghua University, Beijing 100084, China}
\author{Yong-Chun Liu}
\email{ycliu@tsinghua.edu.cn}

\affiliation{State Key Laboratory of Low-Dimensional Quantum Physics, Department
of Physics, Tsinghua University, Beijing 100084, China}
\affiliation{Frontier Science Center for Quantum Information, Beijing 100084, China}
\begin{abstract}
The non-Hermitian skin effect (NHSE), characterized by boundary-localized
eigenstates under open boundary conditions, represents the key feature
of the non-Hermitian lattice systems. Although the non-Bloch band
theory has achieved success in depicting the NHSE in one-dimensional
(1D) systems, its extension to higher dimensions encounters a fundamental
hurdle in the form of the geometry-dependent skin effect (GDSE), where
the energy spectra and the boundary localization of the eigenstates
rely on the lattice geometry. In this work, we establish the non-Bloch
band theory for two-dimensional (2D) GDSE, by introducing a strip
generalized Brillouin zone (SGBZ) framework. Through taking two sequential
1D thermodynamic limits, first along a major axis and then along a
minor axis, we construct geometry-dependent non-Bloch bands, unraveling
that the GDSE originates from the competition between incompatible
SGBZs. Based on our theory, we derive for the first time a crucial
sufficient condition for the GDSE: the non-Bloch dynamical degeneracy
splitting of SGBZ eigenstates, where a continuous set of degenerate
complex momenta breaks down into a discrete set. Moreover, our SGBZ
formulation reveals that the Amoeba spectrum contains the union of
all possible SGBZ spectra, which bridges the gap between the GDSE
and the Amoeba theory. The proposed SGBZ framework offers a universal
roadmap for exploring non-Hermitian effects in 2D lattice systems,
opening up new avenues for the design of novel non-Hermitian materials
and devices with tailored boundary behaviors.
\end{abstract}
\maketitle
\tableofcontents{}

\section{Introduction}

Non-Hermitian physics, which describes the effective dynamics of open
classical and quantum systems, has emerged as a hot topic in recent
decades \cite{ashidaNonHermitianPhysics2020}. Non-Hermiticity not
only opens an avenue towards novel physics \cite{dingNonHermitianTopologyExceptionalpoint2022,el-ganainyNonHermitianPhysicsPT2018},
but also paves the way for new applications \cite{comaronNonHermitianTopologicalEndmode2020,soneExceptionalNonHermitianTopological2020,longhiBulkEdgeCorrespondence2021,wangNonHermitianMorphingTopological2022,zhuAnomalousSingleModeLasing2022,gaoTwoDimensionalReconfigurableNonHermitian2023,weiNonlinearTopologicalLaser2023,budichNonHermitianTopologicalSensors2020,shaoNonHermitianMoireValley2024,yanAdvancesApplicationsNonHermitian2023}.
Among all non-Hermitian systems, non-Hermitian periodic lattices have
attracted intense interest because of their strong connection to the
topological band theory \cite{hasanColloquiumTopologicalInsulators2010,bansilColloquiumTopologicalBand2016}.
In non-Hermitian lattices, the 10-fold Altland-Zirnbauer symmetry
classes \cite{altlandNonstandardSymmetryClasses1997} are generalized
to 38-fold classes \cite{kawabataSymmetryTopologyNonHermitian2019,ashidaNonHermitianPhysics2020,bergholtzExceptionalTopologyNonHermitian2021,dingNonHermitianTopologyExceptionalpoint2022},
providing exotic new topology phases without Hermitian counterparts.

Despite the intriguing properties exhibited by non-Hermitian lattices,
the language of the Bloch band theory, which describes the lattices
under the OBC by the energy bands under the periodic boundary condition
(PBC), breaks down in non-Hermitian lattices. Non-Hermitian lattices
can exhibit the NHSE \cite{leeAnomalousEdgeState2016,yaoEdgeStatesTopological2018},
where the spectrum of the system under the OBC differs greatly from
the spectrum under the periodic boundary condition PBC, and most (proportional
to the size of the system) of the eigenstates under OBC are localized
at the boundaries. Furthermore, when NHSE occurs, topological numbers
defined on the Brillouin zone (BZ) are not consistent with the number
of the topological edge states. Until now, NHSE has been observed
in classical systems such as photonic crystals \cite{weidemannTopologicalFunnelingLight2020,wangGeneratingArbitraryTopological2021,liuComplexSkinModes2022,liuLocalizationChiralEdge2024},
phononic metamaterials \cite{zhangObservationHigherorderNonHermitian2021,guTransientNonHermitianSkin2022,wangExtendedStateLocalized2022}
and topolectrical circuits \cite{helbigGeneralizedBulkBoundary2020,hofmannReciprocalSkinEffect2020,zouObservationHybridHigherorder2021,shangExperimentalIdentificationSecondOrder2022,liuExperimentalObservationNonHermitian2023,wuEvidencingNonBlochDynamics2023,zhangElectricalCircuitRealization2023,zhuHigherRankChirality2023},
as well as in quantum platforms like superconductors \cite{chuBroadColossalEdge2023},
solid-state systems \cite{zhangObservationNonHermitianTopology2021,ochkanNonHermitianTopologyMultiterminal2024}
and photonic quantum-walk platforms \cite{xiaoNonHermitianBulkBoundary2020,wangDetectingNonBlochTopological2021,xiaoObservationNonBlochParityTime2021,linObservationNonHermitianTopological2022,linTopologicalPhaseTransitions2022}.

After the discovery of the NHSE, many theoretical works are proposed
to understand this phenomenon, such as the theory of spectral topology
\cite{leykamEdgeModesDegeneracies2017,gongTopologicalPhasesNonHermitian2018,kunstBiorthogonalBulkBoundaryCorrespondence2018,shenTopologicalBandTheory2018,xiongWhyDoesBulk2018,yinGeometricalMeaningWinding2018,kawabataSymmetryTopologyNonHermitian2019}
and the transfer-matrix approach \cite{kunstNonHermitianSystemsTopology2019}.
However, these theories cannot give a quantitative description of
the thermodynamic limit of OBC spectra and eigenstates just like BZ
in Hermitian lattices. In 1D lattices, this problem is solved by the
non-Bloch band theory \cite{yaoEdgeStatesTopological2018,yokomizoNonBlochBandTheory2019,yangNonHermitianBulkBoundaryCorrespondence2020}.
The main idea of the non-Bloch band theory is to generalize the Bloch
wavevectors to complex numbers, so that the generalized Bloch waves
with imaginary wavevectors (named ``non-Bloch waves'') exhibit exponential
decay in real space. By pushing the OBC to the thermodynamic limit,
the phase factor $\beta=\exp\left(\mu+i\theta\right)\in\mathbb{C}$
of the complex wavevector $k=\theta-i\mu$ is confined onto a 1D closed
loop dubbed ``generalized Brillouin zone'' (GBZ) \cite{yaoEdgeStatesTopological2018}.
Similar to the BZ in Hermitian lattices, the GBZ gives the non-Bloch
energy bands for non-Hermitian lattices, which describes the eigensystem
under OBC.

With the non-Bloch band theory, the anomalous phenomena associated
with the NHSE, i.e., the deviation of the OBC spectrum from the PBC
spectrum, the boundary localization of the OBC eigenstates and the
breakdown of bulk-boundary correspondence, can be quantitatively described
or remedied in 1D lattices. For the spectra, the thermodynamic limit
of the OBC spectrum tends to the spectrum of the GBZ. For the eigenstates,
the exponential localization of eigenstates is described by the imaginary
part of the complex wavevector, or equivalently the modulus of the
phase factor $\beta$ on the GBZ. For the bulk-boundary correspondence,
the existence of edge states is consistent with the topological numbers
calculated on the GBZ. The recovered bulk-boundary correspondence,
dubbed ``non-Bloch bulk boundary correspondence'', is manifested
in various theoretical \cite{okumaTopologicalPhaseTransition2019,zhuPhotonicNonHermitianSkin2020,bartlettIlluminatingBulkboundaryCorrespondence2021,longhiNonHermitianTopologicalPhase2021}
and experimental \cite{xiaoNonHermitianBulkBoundary2020,wangDetectingNonBlochTopological2021,xiaoObservationNonBlochParityTime2021}
researches. Besides the description of OBC spectrum and the recovery
of bulk-boundary correspondence, non-Bloch theory is also applied
in describing the non-Hermitian dynamical evolution \cite{longhiNonBlochBandCollapseChiral2020,xueSimpleFormulasDirectional2021,liNonBlochQuenchDynamics2021,longhiSelfHealingNonHermitianTopological2022,huGreensFunctionsMultiband2023,chenFormalGreensFunction2024}
and the non-Hermitian band engineering \cite{zhouEngineeringNonHermitianSkin2022,jiangReciprocatingBipolarNonHermitian2023,keFloquetEngineeringNonHermitian2023}.
It is also generalized to new systems such as continuum systems \cite{longhiNonHermitianSkinEffect2021,yokomizoNonHermitianPhysicsLevitated2023,yokomizoNonBlochBandTheory2024},
Floquet systems \cite{zhouDualTopologicalCharacterization2021,liLossinducedFloquetNonHermitian2023},
domain wall systems \cite{dengNonBlochTopologicalInvariants2019}
and disordered systems \cite{okumaNonHermitianSkinEffects2021,liuModifiedGeneralizedBrillouin2023,zhangBulkboundaryCorrespondenceDisordered2023}.

After the establishment of the 1D non-Bloch theory, various attempts
are reported to generalize the non-Bloch theory to 2D or higher dimensional
non-Hermitian systems \cite{yaoNonHermitianChernBands2018,yokomizoNonBlochBandsTwodimensional2023,jiangDimensionalTransmutationNonHermiticity2023,wangAmoebaFormulationNonBloch2024,zhangEdgeTheoryNonHermitian2024,huTopologicalOriginNonHermitian2024}.
The most notable one of them is the Amoeba formulation, which defines
the GBZ in arbitrary dimension with a mathematical concept dubbed
``Amoeba'' \cite{wangAmoebaFormulationNonBloch2024}. However, all
of these works fail to describe a special type of NHSE in 2D or higher-dimensional
systems termed GDSE \cite{zhangUniversalNonHermitianSkin2022}. As
schematically illustrated in Fig. \ref{fig:schematic-of-DGBZ}, in
2D or higher-dimensional non-Hermitian systems, the OBC spectrum and
OBC eigenstates can be dependent on the spatial geometry even in the
thermodynamic limit. When a non-Hermitian lattice is truncated into
two different geometries ($G_{1}$ and $G_{2}$ in Fig. \ref{fig:schematic-of-DGBZ}),
the OBC spectra under the two different geometries are different ($\sigma_{1}$
and $\sigma_{2}$ in Fig. \ref{fig:schematic-of-DGBZ}). Moreover,
the eigenstates can exhibit different exponents of boundary localization
in different geometries. For example, the eigenstate of one geometry
($G_{1}$) can be extended in the bulk, shown as $\psi_{1}$ in the
bottom-left of Fig. \ref{fig:schematic-of-DGBZ}, while the eigenstate
with the same energy but under another geometry ($G_{2}$) is localized
at the boundary, shown as $\psi_{2}$ in the bottom-right of Fig.
\ref{fig:schematic-of-DGBZ}. Until now, the GDSE has been observed
in various systems theoretically \cite{chengTruncationdependent$mathcalPT$Phase2022,fangGeometrydependentSkinEffects2022,wangNonHermitianOpticalAtomic2022,sarkarNonHermiticityInducedExceptional2023,qinGeometrydependentSkinEffect2024}
and experimentally \cite{wanObservationGeometrydependentSkin2023,PhysRevLett.131.207201,zhouObservationGeometrydependentSkin2023}.

\begin{figure}
\centering

\includegraphics{Figures/GDSE-20250402}

\caption{Schematic diagram of GDSE. For the same 2D non-Hermitian lattice,
different geometries (blue and orange regions marked by $G_{1}$ and
$G_{2}$) will lead to different spectra ($\sigma_{1}$ and $\sigma_{2}$
in the top-left and top-right plots) and eigenstates ($\psi_{1}$
and $\psi_{2}$ in the bottom-left and bottom-right plots).}
\label{fig:schematic-of-DGBZ} 
\end{figure}

Apart from the non-Bloch band theories, some researchers focus on
the Bloch wave dynamics rather than the non-Bloch waves in 2D and
higher-dimensional non-Hermitian lattices \cite{zhangDynamicalDegeneracySplitting2023,wanObservationGeometrydependentSkin2023}.
They investigate the reflection of Bloch waves at the boundaries,
and find that the localization of the reflected wave is related to
a phenomenon named ``dynamical degenerate splitting'' (DDS), which
depends on the direction of the incident wave. However, Bloch wave
dynamics cannot give the OBC spectrum and the eigenstates, so that
it cannot determine whether the GDSE occurs in a specific non-Hermitian
lattice.

Until now, it is still an enormous challenge to build a theoretical
description for 2D non-Hermitian lattices that suffices to give a
thorough understanding and a universal condition for the GDSE. To
accomplish this task, we derive the strip generalized Brillouin zone
(SGBZ) formulation by taking two sequential 1D thermodynamic limits,
first along a major axis and then along a minor axis. With theoretical
and numerical analysis, it is verified that the SGBZ describes OBC
eigensystems under a strip geometry. When the system exhibits more
than one inequivalent SGBZs, competitions between SGBZs will occur,
implying the existence of the GDSE. With the tool of the SGBZ, we
successfully derive the condition for the GDSE. By comparing the SGBZs
in different strips, we find that the degeneracy splitting of the
non-Bloch waves from a continuum subset of the SGBZ to a discrete
subset, dubbed “non-Bloch DDS” in this work, is a sufficient condition
for the GDSE. Compared with ref. \cite{zhangDynamicalDegeneracySplitting2023},
this non-Bloch DDS is a generalization of the DDS of Bloch waves,
and it will reduce to the latter when the SGBZ coincides with the
BZ. Furthermore, our theory also explains the contradiction between
the GDSE and the Amoeba theory. We prove that the Amoeba spectrum
can be viewed as a combination of the SGBZ spectra in all possible
geometries of a system.

This article is organized as follows. The exact form of SGBZ is derived
in Sec. \ref{sec:strip-GBZ}, and its relation with the GDSE is investigated
with a specific example in Sec. \ref{sec:OBC}. In Sec. \ref{sec:coord-trans},
the transformations of SGBZs are investigated in general to derive
the condition for the GDSE. The non-Bloch DDS is also discussed. Finally,
in Sec. \ref{sec:Amoeba}, the relation between the SGBZ and the Amoeba
formulation is investigated.

\section{SGBZ formulation for 2D non-Hermitian lattices}

\label{sec:strip-GBZ}

In this section, we will build the geometry-dependent non-Bloch band
theory in 2D lattices with the SGBZ formulation. We will first review
the main idea of the non-Bloch band theory, then establish the SGBZ
formulation.

\subsection{Brief introduction to the non-Bloch band theory and 1D GBZ constraint}

Before deriving the SGBZ, we first go through the non-Bloch band theory
and introduce the 1D GBZ constraint as a preparation. For a general
$D$-dimensional non-interacting periodic lattice, with a basis of
lattice vectors $\left\{ \mathbf{a}_{1},\mathbf{a}_{2},\dots,\mathbf{a}_{D}\right\} $,
the Hamiltonian of the system can be written in a general form, which
reads, 
\begin{equation}
H=\sum_{\mathbf{r}}\sum_{t_{1},t_{2},\dots,t_{D}}\sum_{\mu,\nu=1}^{m}\mathcal{T}_{\mathbf{t},\mu,\nu}c_{\mathbf{r}+\mathbf{t},\mu}^{\dagger}c_{\mathbf{r},\nu},\label{eq:non-interacting-Hamiltonian}
\end{equation}
where $\mathcal{T}_{\mathbf{t},\mu,\nu}\in\mathbb{C}$ is the coupling
coefficient, and $c_{\mathbf{r},\nu}$ is the annihilator at position
$\mathbf{r}$ and sublattice index $\nu$. The subscript $\mathbf{r}=\sum_{j=1}^{D}r_{j}\mathbf{a}_{j}$
is the position vector, and $\mathbf{t}=\sum_{j=1}^{D}t_{j}\mathbf{a}_{j}$
is the coupling vector within the range $|t_{j}|\leq t_{Rj},j=1,2,\dots,D$,
where $t_{Rj},j=1,2,\dots,D$ are positive integers. The coordinates
$r_{j}$ and $t_{j}$ are all integers. The coupling coefficients
satisfy $\mathcal{T}_{\mathbf{t},\mu,\nu}=\mathcal{T}_{-\mathbf{t},\nu,\mu}^{*}$
for Hermitian lattices, and $\mathcal{T}_{\mathbf{t},\mu,\nu}\neq\mathcal{T}_{-\mathbf{t},\nu,\mu}^{*}$
for non-Hermitian lattices.

In Hermitian lattices, according to the Bloch's theorem, the Hamiltonian
$H$ can be diagonalized by the Bloch states due to the discrete translation
symmetries, defined as, 
\begin{equation}
T_{\mathbf{u}}\left|\psi_{\mathbf{k}}\right>=\exp\left(i\sum_{j=1}^{D}u_{j}k_{j}\right)\left|\psi_{\mathbf{k}}\right>,\label{eq:Bloch-states}
\end{equation}
where $T_{\mathbf{u}}$ is the discrete translation operator of the
displacement $\mathbf{u}=u_{1}\mathbf{a}_{1}+u_{2}\mathbf{a}_{2}+\cdots+u_{D}\mathbf{a}_{D}$,
defined as, 
\begin{equation}
T_{\mathbf{u}}c_{\mathbf{r},\nu}T_{\mathbf{u}}^{-1}=c_{\mathbf{r}+\mathbf{u},\nu}.\label{eq:translation-operators}
\end{equation}
Under the basis of the Bloch states, the Hamiltonian $H$ is block
diagonalized to the momentum-space Hamiltonian $h(e^{i\mathbf{k}})$,
defined as, 
\begin{equation}
h_{\mu,\nu}\left(e^{i\mathbf{k}}\right)=\sum_{t_{1},t_{2},\dots,t_{D}}\mathcal{T}_{\mathbf{t},\mu,\nu}\exp\left(-i\sum_{j=1}^{D}t_{j}k_{j}\right),\label{eq:Bloch-Hamiltonian}
\end{equation}
which is an $m\times m$ matrix. Here, $e^{i\mathbf{k}}$ denotes
the vector $(e^{ik_{1}},e^{ik_{2}},\dots,e^{ik_{D}})$. The good quantum
number $\mathbf{k}$ is called the Bloch wavevector.

In non-Hermitian lattices, due to the existence of the NHSE, the Bloch
waves $|\psi_{\mathbf{k}}\rangle$, which are periodically distributed
in the bulk, cannot describe the boundary localized eigenstates under
OBC. Because the discrete translation symmetries still hold, a direct
remedy of this contradiction is to extend the phase factor $\exp(ik_{j})$
to $\beta_{j}\in\mathbb{C}$, and generalize the Bloch waves $|\psi_{\mathbf{k}}\rangle$
to ``non-Bloch waves'' $|\psi_{\boldsymbol{\beta}}\rangle$, defined
as, 
\begin{equation}
T_{\mathbf{u}}\left|\psi_{\boldsymbol{\beta}}\right>=\left(\prod_{j}\beta_{j}^{u_{j}}\right)\left|\psi_{\boldsymbol{\beta}}\right>,\label{eq:non-Bloch-state}
\end{equation}
where the symbol $\boldsymbol{\beta}$ is defined as $\boldsymbol{\beta}\equiv\left(\beta_{1},\beta_{2},\dots,\beta_{D}\right)$
\cite{yaoEdgeStatesTopological2018}. Similarly, the momentum-space
Hamiltonian $h(e^{i\mathbf{k}})$ is extended to the non-Bloch Hamiltonian
$h(\boldsymbol{\beta})$, which is, 
\begin{equation}
h_{\mu\nu}\left(\boldsymbol{\beta}\right)=\sum_{t_{1},t_{2},\dots,t_{D}}\mathcal{T}_{\mathbf{t},\mu,\nu}\prod_{j=1}^{D}\beta_{j}^{-t_{j}}.\label{eq:non-Bloch-Hamiltonian}
\end{equation}

However, the extension of $e^{i\mathbf{k}},\mathbf{k}\in\mathbb{R}^{D}$
to $\boldsymbol{\beta}\in\mathbb{C}^{D}$ makes the dimensionality
of the momentum space twice as many as the real space, which contradicts
with the numerical results of the OBC spectra of non-Hermitian lattices.
For instance, in 1D non-Hermitian lattices, the OBC spectra are joints
of 1D curves on the complex plane, which contradicts with the non-Bloch
band given by $\beta\in\mathbb{C}$, i.e., $\left\{ E\in\mathbb{C}\mid\det[E-h(\beta)]=0,\beta\in\mathbb{C}\right\} $.
Therefore, additional constraints on $\boldsymbol{\beta}$ are required
to get the physical solutions. In 1D lattices, the constraints are
given by Refs. \cite{yaoEdgeStatesTopological2018,yokomizoNonBlochBandTheory2019,yangNonHermitianBulkBoundaryCorrespondence2020}.
By pushing the OBC to the thermodynamic limit, a real-valued constraint
equation is derived. Assuming that $-M$ and $N$ are the lowest and
highest degrees of $\beta$ in the ChP $f\left(E,\beta\right)\equiv\det[E-h(\beta)]$,
the GBZ constraint reads, 
\begin{equation}
\left|\beta^{(M)}\left(E\right)\right|=\left|\beta^{(M+1)}\left(E\right)\right|,\label{eq:definition-1D-GBZ}
\end{equation}
where $\beta^{(j)}(E)$ is the $j$-th solution of $f(E,\beta)=0$
ordered by $|\beta^{(i)}(E)|\leq|\beta^{(j)}(E)|,\forall i<j$, taking
$E$ as a parameter. When $E$ traverses $\mathbb{C}$, the set of
$\beta^{(M)}(E)$ and $\beta^{(M+1)}(E)$ satisfying Eq. (\ref{eq:definition-1D-GBZ})
is defined as the GBZ. Noted that the GBZ derived from Eq. (\ref{eq:definition-1D-GBZ})
is uniquely determined by the non-Bloch ChP $f(E,\beta)$. Hence,
the GBZ is independent of the size of the system.

Non-Bloch band theory provides a systematic description of the NHSE
in 1D non-Hermitian systems. First, the exponents of the boundary
localization are described by the modulus of $\beta$, because the
profile of the non-Bloch wave is $\left|\langle x|\psi_{\beta}\rangle\right|\sim\left|\beta\right|^{x}$.
When $\left|\beta\right|<1$, the non-Bloch wave is localized at the
left boundary, and vice versa. Second, the thermodynamic limit of
the OBC spectrum is approximated by the GBZ spectrum, i.e., the complex
number $E$ satisfying Eq. (\ref{eq:definition-1D-GBZ}). Third, when
computing the topological numbers on the GBZ instead of the BZ, the
topological numbers are consistent with the number of boundary states.
All the three observations indicate that GBZ is a good substitution
of BZ when a 1D lattice exhibits the NHSE.

However, for 2D and higher-dimensional non-Hermitian lattices, there
is no direct generalization of the 1D non-Bloch band theory. Taking
2D lattices as an example, for the ChP $f(E,\beta_{1},\beta_{2})=\det[E-h(\beta_{1},\beta_{2})]$,
the solutions of the eigenvalue equation $f\left(E,\beta_{1},\beta_{2}\right)=0$
form a continuum on the complex plane rather than a finite set of
points, making $M$ in Eq. (\ref{eq:definition-1D-GBZ}) undefined.
The existence of the GDSE is also an evidence of this obstacle. In
1D non-Hermitian systems, the ``shape'' of a lattice in the thermodynamic
limit is unique, that is, a 1D open chain extended along both sides.
However, in 2D or higher-dimensional systems, the geometry may tend
to infinity in different ways, which cannot be described in the current
version of non-Bloch band theory.

\subsection{Derivation of the SGBZ formulation}

\begin{figure*}
\centering

\includegraphics{Figures/main-idea-20250429}

\caption{Schematics of the SGBZ. (a) Sketch of a 2D lattice. (b) Schematic
diagram of the strip geometry, where $\mathbf{a}_{1}$, $\mathbf{a}_{2}$
denote the major and minor axes, respectively, and the cyan region
denotes the supercell. (c) The hybrid real-momentum space formed by
the 2D complex plane of $\beta_{1}$ and the 1D real space of $r_{2}$.
(d) The QMGBZ, defined as the 1D GBZ of the strip geometry. (e) The
PMGBZ, defined as the parametric GBZ along the minor axis with parameter
$\beta_{1}$. (f) The SGBZ, defined as the limit of the QMGBZ when
the width tends to infinity. The solid arrows indicate direct mathematical
relations between concepts, while the dashed arrows represent indirect
logical relations.}
\label{fig:main-idea} 
\end{figure*}

Owing to the success of 1D non-Bloch band theory, we establish the
2D non-Bloch band theory by imposing two successive 1D GBZ constraints
to $(\beta_{1},\beta_{2})\in\mathbb{C}^{2}$. The main idea is illustrated
in Fig. \ref{fig:main-idea}. For a 2D non-Hermitian lattice {[}Fig.
\ref{fig:main-idea}(a){]}, a strip geometry can be defined by specifying
a basis $\{\mathbf{a}_{1},\mathbf{a}_{2}\}$, and selecting one of
the basis vectors ($\mathbf{a}_{1}$) as the major axis, and the other
($\mathbf{a}_{2}$) as the minor axis, shown as Fig. \ref{fig:main-idea}(b).
The strip geometry is extended along the major axis, and truncated
in the minor axis. When the width of the strip is finite, the 2D lattice
confined in a strip geometry can be viewed as a 1D periodic lattice
with the discrete translation symmetry in the major axis, whose periodic
unit is the supercell marked by the light blue region in Fig. \ref{fig:main-idea}(b).

According to the 1D non-Bloch band theory, when the length in the
major axis tends to infinity, the OBC eigensystem of the strip geometry
tends to the eigensystem of its 1D GBZ, dubbed the quasi-1D major-axis
GBZ (QMGBZ). Supposing the real-space Hamiltonian has the form of
Eq. (\ref{eq:non-interacting-Hamiltonian}), by Fourier transformation
in the major axis, we get a hybrid real-momentum Hamiltonian, which
reads 
\begin{align}
\mathcal{H}\left(\beta_{1}\right)= & \sum_{t_{1},t_{2}}\sum_{r_{2}=1}^{L_{2}}\sum_{\mu,\nu=1}^{m}\mathcal{T}_{t_{1},t_{2},\mu,\nu}\beta_{1}^{-t_{1}}c_{\beta_{1},r_{2}+t_{2};\mu}^{\dagger}c_{\beta_{1},r_{2};\nu},\label{eq:hybrid-non-Bloch-Hamiltonian}
\end{align}
where $L_{2}$ is the width, and $c_{\beta_{1},r_{2};\mu}$ is the
annihilator of the non-Bloch wave with complex momentum of $\exp(ik_{1})\equiv\beta_{1}\in\mathbb{C}$,
the coordinate $r_{2}$ and the sublattice index $\mu$. As illustrated
in Fig. \ref{fig:main-idea}(c), $\beta_{1}$ and $r_{2}$ form a
3D hybrid real-momentum space. To get the QMGBZ {[}Fig. \ref{fig:main-idea}(d){]},
the following 1D GBZ constraint, 
\begin{equation}
|\beta_{1}^{(M)}(E)|=|\beta_{1}^{(M+1)}(E)|,\label{eq:quasi-1D-GBZ-constraint}
\end{equation}
is imposed on the 3D hybrid space, where $\beta_{1}^{(j)}(E)$ is
the $j$-th solution of $F(E,\beta_{1})\equiv\det[E-\mathcal{H}(\beta_{1})]=0$
ordered by $|\beta_{1}^{(i)}(E)|\leq|\beta_{1}^{(j)}(E)|,\forall i<j$,
and $-M,N$ are the lowest and highest degrees of $\beta_{1}$ in
the ChP $F(E,\beta_{1})$, respectively.

Next, to get the SGBZ, we push the width of the strip to infinity
and define the SGBZ as the limit of the QMGBZ. Because QMGBZ is defined
as the points in the hybrid real-momentum space restricted by the
1D GBZ constraint, SGBZ can be derived by the limit of both the hybrid
real-momentum space and the 1D GBZ constraint on $\beta_{1}$ when
the width tends to infinity.

For the hybrid real-momentum space, taking $\beta_{1}$ as a parameter,
the hybrid Hamiltonian $\mathcal{H}(\beta_{1})$ can be viewed as
a parametric Hamiltonian of a 1D open chain along the minor axis.
When the width tends to infinity, the eigensystem of $\mathcal{H}(\beta_{1})$
is consistent with the parametric 1D GBZ along the minor axis, called
the parametric minor-axis GBZ (PMGBZ). As illustrated in Fig. \ref{fig:main-idea}(e),
each $\beta_{1}\in\mathbb{C}$ specifies a 1D GBZ of $\beta_{2}$.
When $\beta_{1}$ traverses $\mathbb{C}$, all these 1D GBZs form
a 3D space.

For the constraint, because the degree $M$ in Eq. (\ref{eq:quasi-1D-GBZ-constraint})
diverges when the width tends to infinity, we use the winding number
formulation for the 1D GBZ constraint \cite{zhangCorrespondenceWindingNumbers2020},
which is based on the quasi-1D winding number defined as, 
\begin{equation}
W\left(E,r\right)=\frac{1}{2\pi i}\oint_{C_{r}}d\beta_{1}\frac{d\ln F\left(E,\beta_{1}\right)}{d\beta_{1}},\label{eq:definition-winding-number-quasi-1D}
\end{equation}
where $C_{r}\equiv\left\{ \beta_{1}\in\mathbb{C}\mid\left|\beta_{1}\right|=r\right\} $
is the circle centered at $0$ with radius $r$. When the width tends
to infinity, we prove that $W(E,r)/L_{2}$ converges to a quantity
$W_{\text{strip}}(E,r)$, named the ``strip winding number'', which
is uniquely determined by the 2D non-Bloch Hamiltonian $h(\beta_{1},\beta_{2})$
defined as Eq. (\ref{eq:non-Bloch-Hamiltonian}). Substituting $W_{\text{strip}}(E,r)$
into $W(E,r)$ in the constraint for QMGBZ, we get the constraint
for the SGBZ.

In the following two parts, we will discuss the details about the
winding number formulation of the 1D GBZ constraint and the strip
winding number $W_{\text{strip}}(E,r)$, respectively.

\subsubsection{Winding number formulation and its thermodynamic limit}

\label{subsec:Winding-number-formulation}

For the quasi-1D winding number of Eq. (\ref{eq:definition-winding-number-quasi-1D}),
when $F(E,\beta_{1})$ does not vanish on $C_{r}$, according to the
Cauchy's argument principle, the quasi-1D winding number is related
to the $\beta_{1}$-solutions of $F(E,\beta_{1})=0$ by, 
\begin{equation}
W\left(E,r\right)=N_{\text{zeros}}-M,\label{eq:winding-number-and-zeros-1}
\end{equation}
where $N_{\text{zeros}}$ is the number of $\beta_{1}$-solutions
enclosed by $C_{r}$. Therefore, if the constraint for QMGBZ is not
satisfied for some reference energy $E$, that is, $|\beta_{1}^{(M)}(E)|<|\beta_{1}^{(M+1)}(E)|$,
the quasi-1D winding number $W(E,r)$ vanishes when $|\beta_{1}^{(M)}(E)|<r<|\beta_{1}^{(M+1)}(E)|$.
Based on this observation, the 1D GBZ constraint is equivalent to
the following form: for some reference energy $E$ and radius $r$,
if there exist $r_{<}\in(r-\epsilon,r)$ and $r_{>}\in(r,r+\epsilon)$
for arbitrary $\epsilon>0$, such that, 
\begin{equation}
\begin{cases}
W(E,r_{<})<0,\\
W(E,r_{>})>0,
\end{cases}\label{eq:quasi-1D-winding-number-constraint}
\end{equation}
then the $\beta_{1}$-solutions of $F(E,\beta_{1})=0$ satisfying
$|\beta_{1}|=r$ are in the QMGBZ, and vice versa. The geometry picture
of the winding number formulation is shown in Fig. \ref{fig:q1D-thermo}(a).
According to Eq. (\ref{eq:winding-number-and-zeros-1}) and Eq. (\ref{eq:quasi-1D-winding-number-constraint}),
the winding loop $C_{r_{<}}$ encloses no more than $M-1$ zeros,
and the loop $C_{r_{>}}$ encloses no less than $M+1$ zeros. Therefore,
$|\beta_{1}^{(M)}(E)|$ and $|\beta_{1}^{(M+1)}(E)|$ belong to $(r_{<},r_{>})\subset(r-\epsilon,r+\epsilon)$
for arbitrary $\epsilon$. When $\epsilon$ reduces to $0$, we have
$|\beta_{1}^{(M)}(E)|=|\beta_{1}^{(M+1)}(E)|=r$, which is equivalent
to Eq. (\ref{eq:quasi-1D-GBZ-constraint}).

\begin{figure}
\centering

\includegraphics{Figures/thermodynamic-quasi-1D-20250505}

\caption{The winding number formulation of the 1D GBZ constraint, and the thermodynamic
distribution of the $\beta_{1}$-zeros. (a, b) Illustration of the
1D GBZ constraint for the QMGBZ with (a) finite width and (b) infinite
width. The blue points or curves denote the $\beta_{1}$-zeros of
$F(E,\beta_{1})$, and the green dots denote the QMGBZ solutions.
At thermodynamic limit, the $\beta_{1}$-zeros tend to $\text{PMGBZ}(E)$.
The red and purple dotted circles denote the winding loops with radii
$r_{<}$ and $r_{>}$, respectively. (c) Distribution of $\beta_{1}$-zeros
on a segment of $\text{PMGBZ}(E)$, where $\theta_{2}\equiv\text{Arg}(\beta_{2})$.
The relative phases at the both endpoints of the segment are marked
by the red arrows. (d) Distribution of the $\beta_{1}$-zeros between
two circular winding loops. The black dashed circles represent the
winding loops $C_{r}$ and $C_{r^{\prime}}$, respectively. The blue
curves denote the segments of $\text{PMGBZ}(E)$ between the two winding
loops, and the purple curves denote its projection on $\beta_{1}$-plane.
The red arrows mark the relative phase.}
\label{fig:q1D-thermo} 
\end{figure}

To obtain the thermodynamic limit of $W(E,r)$, the thermodynamic
distribution of the zeros $\beta_{1}^{(j)}(E)$ is required. By definition,
the ChP $F(E,\beta_{1})$ can be expanded as, 
\begin{equation}
F(E,\beta_{1})\equiv\det\left[E-\mathcal{H}(\beta_{1})\right]=\prod_{j=1}^{mL_{2}}\left(E-E_{j}(\beta_{1})\right),\label{eq:QMGBZ-determinant}
\end{equation}
where $E_{j}(\beta_{1}),j=1,2,\dots,mL_{2}$ are the eigenvalues of
$\mathcal{H}(\beta_{1})$. Therefore, for a given reference energy
$E$, the $\beta_{1}$-zeros of $F(E,\beta_{1})$ are the values of
$\beta_{1}$ where the spectrum of $\mathcal{H}(\beta_{1})$ contains
$E$. As discussed above, when the width tends to infinity, the spectrum
of $\mathcal{H}(\beta_{1})$ tends to the spectrum of the 1D GBZ in
the minor axis, called the PMGBZ. According to the 1D non-Bloch band
theory, the PMGBZ are the set of the points $\left(\beta_{1},\beta_{2}^{(M_{2})}(E,\beta_{1})\right)$
and $\left(\beta_{1},\beta_{2}^{(M_{2}+1)}(E,\beta_{1})\right)$ that
satisfy the 1D GBZ constraint in the minor axis, 
\begin{equation}
\left|\beta_{2}^{(M_{2})}(E,\beta_{1})\right|=\left|\beta_{2}^{(M_{2}+1)}(E,\beta_{1})\right|,\label{eq:PMGBZ-constraint}
\end{equation}
where $\beta_{2}^{(j)}(E,\beta_{1})$ is the solution of $f(E,\beta_{1},\beta_{2})\equiv\det[E-h(\beta_{1},\beta_{2})]=0$
ordered by $|\beta_{2}^{(i)}(E,\beta_{1})|\leq|\beta_{2}^{(j)}(E,\beta_{1})|,\forall i<j$,
and $-M_{2}$ is the lowest degree of $\beta_{2}$ in $f(E,\beta_{1},\beta_{2})$.
The thermodynamic distribution of the zeros $\beta_{1}^{(j)}(E)$
is the PMGBZ points with eigenenergy $E$, i.e., 
\begin{equation}
\text{PMGBZ}\left(E\right)=\left\{ \left(\beta_{1},\beta_{2}\right)\in\text{PMGBZ}\mid f\left(E,\beta_{1},\beta_{2}\right)=0\right\} .\label{eq:PMGBZ-of-E}
\end{equation}
As sketched in Fig. \ref{fig:q1D-thermo}(b), the set $\text{PMGBZ}(E)$
is locally 1D for a fixed reference energy $E$, because it is a subset
of a 4D space $(\beta_{1},\beta_{2})\in\mathbb{C}^{2}$ restricted
by a complex-valued constraint {[}eigenvalue equation $f(E,\beta_{1},\beta_{2})=0${]}
and a real-valued constraint {[}Eq. (\ref{eq:PMGBZ-constraint}){]}.

According to Eq. (\ref{eq:PMGBZ-constraint}), points in PMGBZ are
in pairs. For a PMGBZ pair $(\beta_{1},\beta_{2}^{(\mathrm{a})})$
and $(\beta_{1},\beta_{2}^{(\mathrm{b})})$ in $\text{PMGBZ}(E)$,
we prove that the relative phase $\phi\equiv\text{Arg}(\beta_{2}^{(\mathrm{b})}/\beta_{2}^{(\mathrm{a})})$
is related to the density of $\beta_{1}$-zeros of $F(E,\beta_{1})$
(see Sec. \ref{app:DOS-on-1D-GBZ} of the Supplementary Materials
(SM) \cite{supp-2SK}). As illustrated in Fig. \ref{fig:q1D-thermo}(c),
when $L_{2}$ is large enough, the relation between the $\beta_{1}$-zeros
of $F(E,\beta_{1})$ and the relative phase reads, 
\begin{equation}
N_{\text{q1D}}=\frac{L_{2}}{2\pi}|\phi-\phi^{\prime}|+O(1),
\end{equation}
where $N_{\text{q1D}}$ is the number of the $\beta_{1}$-zeros located
on the projection of the PMGBZ segment on the $\beta_{1}$-plane,
and $\phi$, $\phi^{\prime}$ are the relative phases at the endpoints
of the pair of segments \footnote{When $r^{\prime}$ is close enough to $r$, this selection of base
points is always valid except for the pairs with zero relative phase,
which consist a set with Lebesgue measure zero on $\text{mGBZ}(E)$.}. Therefore, we obtain the thermodynamic distribution of $W(E,r)$.
For two nearby circles with radii $r^{\prime}$ and $r$ ($r^{\prime}>r$),
if the two circles intersect with $k$ pairs of PMGBZ segments, as
sketched in Fig. \ref{fig:q1D-thermo}(d), the increment of $W(E,r)$
reads, 
\begin{equation}
W(E,r^{\prime})-W(E,r)=\frac{L_{2}}{2\pi}\sum_{j=1}^{k}\left|\phi_{j}-\phi_{j}^{\prime}\right|+O\left(1\right),\label{eq:increment-q1D-winding-number}
\end{equation}
where $\phi_{j}$ and $\phi_{j}^{\prime}$ denote the relative phases
of the $j$-th pair of segments at $C_{r}$ and $C_{r^{\prime}}$,
respectively.

\subsubsection{Strip winding number and definition of SGBZ}

\label{subsec:Strip-winding-number}

Inspired by Eq. (\ref{eq:increment-q1D-winding-number}), we consider
the limit of $W(E,r)/L_{2}$ when $L_{2}\rightarrow\infty$. We find
that $W(E,r)/L_{2}$ converges to a quantity $W_{\text{strip}}(E,r)$
named the ``strip winding number'', which is uniquely determined
by the non-Bloch Hamiltonian $h(\beta_{1},\beta_{2})$. We will first
give the definition of $W_{\text{strip}}(E,r)$, and then check its
relation to $W(E,r)$.

\begin{figure}
\centering

\includegraphics{Figures/total-winding-number-20250414}

\caption{Sketches for the definition and the properties of the strip winding
number. (a) Definition of the minor radius function $\rho(\theta_{1})$,
where the blue solid line and the cyan dashed line denote the $M_{2}$-th
and $(M_{2}+1)$-th $\beta_{2}$-zeros of $f\left(E,re^{i\theta_{1}},\beta_{2}\right)$.
(b) Illustration of the winding number $w(\theta_{2})$ and PMGBZ
points on the base manifold $X(E,r)$. The cyan dashed lines and blue
solid lines marked by $(M_{2})$ and $(M_{2}+1)$ denote the $M_{2}$-th
and $(M_{2}+1)$-th $\beta_{2}$-zeros, the orange curves denote the
winding loops of $w(\theta_{2})$, the red and green points denote
the PMGBZ points with topological charges $+1$ and $-1$, respectively,
and the loop $\ell_{0}$ denote a horizontal loop encircling a curve
of $\beta_{2}$-zeros. (c) Relation of $w(\theta_{2})$ and PMGBZ
pairs on the unfolded manifold, where the red and green points denote
the PMGBZ pairs, and the orange arrowed line denote a special winding
loop. The values of $w(\theta_{2})$ in different colored stripes
are marked in the bottom of each region. (d) Change of PMGBZ pairs
when $|\beta_{1}|$ increases. When $|\beta_{1}|$ increases from
$r$ to $r^{\prime}>r$, the distance from \textquotedblleft$+$\textquotedblright{}
to \textquotedblleft$-$\textquotedblright{} increases, and the distance
from \textquotedblleft$-$\textquotedblright{} to \textquotedblleft$+$\textquotedblright{}
decreases.}
\label{fig:strip-winding} 
\end{figure}

For a given reference energy $E$ and a radius $r$, we consider the
$\beta_{2}$-zeros of 2D ChP $f(E,\beta_{1},\beta_{2})$, with $\beta_{1}=re^{i\theta_{1}}$
to be a parameter. As illustrated in Fig. \ref{fig:strip-winding}(a),
we pick a $2\pi$-periodic function $\rho(\theta_{1})$ sandwiched
between the moduli of the $M_{2}$-th and $(M_{2}+1)$-th $\beta_{2}$-zeros
of $f(E,\beta_{1},\beta_{2})$, that is, 
\begin{equation}
\left|\beta_{2}^{(M_{2})}\left(E,re^{i\theta_{1}}\right)\right|\leq\rho\left(\theta_{1}\right)\leq\left|\beta_{2}^{(M_{2}+1)}\left(E,re^{i\theta_{1}}\right)\right|,
\end{equation}
where ``$=$'' holds if and only if $|\beta_{2}^{(M_{2})}|=|\beta_{2}^{(M_{2}+1)}|$.
To define $W_{\text{strip}}(E,r)$, as illustrated in Fig. \ref{fig:strip-winding}(b),
a 2D manifold $X(E,r)$ is constructed as, 
\begin{equation}
X\left(E,r\right)\equiv\left\{ \left(re^{i\theta_{1}},\rho\left(\theta_{1}\right)e^{i\theta_{2}}\right)\mid\theta_{1},\theta_{2}\in\left[-\pi,\pi\right]\right\} .\label{eq:toric-base-manifold}
\end{equation}
Owing to the periodicity of $\rho(\theta_{1})$, the manifold of $X(E,r)$
has the topology of a 2D torus. According to the definition of $\rho(\theta_{1})$,
the zeros $\beta_{2}^{(M_{2})}$ (cyan dashed curves) are located
inside the torus $X(E,r)$, and $\beta_{2}^{(M_{2}+1)}$ (blue solid
curves) are outside of $X(E,r)$. The intersections between the curves
of the $\beta_{2}$-zeros and the manifold $X(E,r)$ are PMGBZ pairs,
shown as the red and green dots in Fig. \ref{fig:strip-winding}(b).

With the above preparations, the strip winding number is defined as
the integration, 
\begin{equation}
W_{\text{strip}}\left(E,r\right)=\frac{1}{2\pi}\int_{-\pi}^{\pi}w\left(\theta_{2}\right)d\theta_{2},\label{eq:definition-W-strip}
\end{equation}
where $w(\theta_{2})$ is the winding number defined as, 
\begin{equation}
w\left(\theta_{2}\right)=\int_{0}^{2\pi}\frac{d\theta_{1}}{2\pi i}\frac{\partial\ln\left[f\left(E,re^{i\theta_{1}},\rho\left(\theta_{1}\right)e^{i\theta_{2}}\right)\right]}{\partial\theta_{1}}.\label{eq:loop-winding-number}
\end{equation}
From the geometric perspective, $w(\theta_{2})$ defined in Eq. (\ref{eq:loop-winding-number})
is the winding number of ChP $f(E,re^{i\theta_{1}},\rho(\theta_{1})e^{i\theta_{2}})$
when $\theta_{1}$ increases by $2\pi$ and $\theta_{2}$ remains
constant, as illustrated by the orange curves in Fig. \ref{fig:strip-winding}(b).
Because of the topological invariance, $w(\theta_{2})$ keeps constant
under the homotopic deformation of the winding loop that keeps the
ChP non-vanishing. On the base manifold $X(E,r)$, the ChP vanishes
and only vanishes at the PMGBZ points, so that $w(\theta_{2})$ changes
if and only if the winding loop passes the PMGBZ points.

When the winding loop of $w(\theta_{2})$ moves across an PMGBZ point
along the positive direction of $\theta_{2}$, the increment of $w(\theta_{2})$
equals the winding number around an infinitesimal closed loop encircling
the PMGBZ point. The direction of the infinitesimal loop is right-handed
to the outward normal vector of $X(E,r)$, shown as the red and green
loops around the PMGBZ points in Fig. \ref{fig:strip-winding}(b).
A detailed discussion is available in Sec. \ref{app:winding-general-loop}
of the SM \cite{supp-2SK}. Because the winding number of the infinitesimal
loop is determined by the local properties of a PMGBZ point, it can
be viewed as a topological charge attached to the PMGBZ point. Furthermore,
the topological charge is determined by the direction in which the
curve of the $\beta_{2}$-zeros passes $X(E,r)$. As illustrated in
Fig. \ref{fig:strip-winding}(b), consider the horizontal loop $\ell_{0}$
encircling one of the $\beta_{2}$-zero points $\beta_{2}^{(M_{2})}$
or $\beta_{2}^{(M_{2}+1)}$, parameterized as, 
\begin{equation}
\begin{cases}
\theta_{1}(t)=\theta_{1,\ell_{0}}\\
\beta_{2}(t)=\beta_{2}^{(j)}\left(E,re^{i\theta_{1,\ell_{0}}}\right)+\epsilon e^{it}
\end{cases}t\in\left[0,2\pi\right],\label{eq:parameterized-sigma-0}
\end{equation}
where $\theta_{1,\ell_{0}}$ is an arbitrary constant, and $j=M_{2}$
or $M_{2}+1$ is the index of the $\beta_{2}$-zero. When $\epsilon$
is small enough, the ChP reads, 
\begin{equation}
f\left(E,re^{i\theta_{1}(t)},\beta_{2}\left(t\right)\right)=\epsilon\frac{\partial f}{\partial\beta_{2}}e^{it}+O\left(\epsilon^{2}\right),\label{eq:expansion-ChP-sigma-0}
\end{equation}
where the derivative $\partial f/\partial\beta_{2}$ is evaluated
at the zero $\left(E,re^{i\theta_{1,\ell_{0}}},\beta_{2}^{(j)}(E,re^{i\theta_{1,\ell_{0}}})\right)$.
When $\partial f/\partial\beta_{2}\neq0$, the winding number around
$\ell_{0}$ reads, 
\begin{align}
w_{\text{loop}}\left(\ell_{0}\right) & \equiv\int_{0}^{2\pi}\frac{dt}{2\pi i}\frac{\partial\ln\left[f\left(E,re^{i\theta_{1}(t)},\beta_{2}\left(t\right)\right)\right]}{\partial t},\nonumber \\
 & =1.\label{eq:winding-number-sigma-0}
\end{align}
Therefore, all loops homotopic to $\ell_{0}$ take a winding number
of $+1$, and all loops homotopic to the reverse of $\ell_{0}$ take
a winding number of $-1$. As illustrated in Fig. \ref{fig:strip-winding}(b),
when $\theta_{1}$ increases, if the $\beta_{2}$-solution moves outward
$X(E,r)$ at a PMGBZ point, the winding loop around the PMGBZ point
is homotopic to $\ell_{0}$, so that the topological charge of the
PMGBZ point is $+1$, and vice versa. We also conclude that the topological
charges of a pair of PMGBZ points are opposite to each other, according
to the definition of the radius function $\rho(\theta_{1})$.

Due to the topological invariance of $w(\theta_{2})$, the expression
of $W_{\text{strip}}(E,r)$ can be simplified. In Fig. \ref{fig:strip-winding}(c),
the PMGBZ pairs are shown in the unfolded view of $X(E,r)$. As illustrated
by the colored stripes in Fig. \ref{fig:strip-winding}(c), the PMGBZ
points split $X(E,r)$ into different slices parallel to the $\theta_{1}$-axis.
In each slice, $w(\theta_{2})$ is constant. We pick an arbitrary
loop with $\theta_{2}=\theta_{2,0}$, shown as the orange arrowed
line in Fig. \ref{fig:strip-winding}(c), and calculate its winding
number $w(\theta_{2,0})=w_{0}$. Then, the values of $w(\theta_{2})$
in other slices are determined according to the topological charges
of the PMGBZ points, marked at the bottom of each slice in Fig. \ref{fig:strip-winding}(c).
As an example, for the case illustrated in Fig. \ref{fig:strip-winding}(c),
the strip winding number reads $W_{\text{strip}}(E,r)=w_{0}-\phi_{1}/2\pi+\phi_{2}/2\pi$,
where $\phi_{1}$ and $\phi_{2}$ are the relative phases of the two
PMGBZ pairs, respectively.

In general, if there exist $k$ pairs of PMGBZ points on the manifold
$X(E,r)$, we first pick an arbitrary winding loop with $\theta_{2}=\theta_{2,0}$,
then calculate its winding number $w(\theta_{2,0})=w_{0}$. Each pair
of PMGBZ points cut the circle parallel to the $\theta_{2}$-axis
{[}dotted lines in Fig. \ref{fig:strip-winding}(c). Noted that the
horizontal lines in Fig. \ref{fig:strip-winding}(c) are circles because
of the periodicity of $\theta_{2}${]} into two different arcs. Similar
to Fig. \ref{fig:strip-winding}(c), for $j$-th PMGBZ point, we pick
the arc that does not intersect with the line $\theta_{2}=\theta_{2,0}$,
and denote the relative phase on the selected arc as $\phi_{j}$.
Then, the winding number of the loop that intersects the selected
arc is larger or smaller than the loop that does not intersect the
selected arc by $1$, where the sign of the increment is determined
by the topological charges. Therefore, the $j$-th PMGBZ pair contributes
to $W_{\text{strip}}(E,r)$ by $\pm\phi_{j}/2\pi$, and the general
expression of the strip winding number reads, 
\begin{equation}
W_{\text{strip}}\left(E,r\right)=w_{0}+\sum_{j=1}^{k}\left(-1\right)^{\tau_{j}}\frac{\phi_{j}}{2\pi},\label{eq:W-strip-general}
\end{equation}
where $\tau_{j}=0,1$ depends on the topological charges of the $j$-th
PMGBZ pair. When the start point of the selected arc has a positive
charge, $\tau_{j}=0$, and vice versa.

With the simplified relation of Eq. (\ref{eq:W-strip-general}), we
will verify the relation between $W(E,r)$ and $W_{\text{strip}}(E,r)$.
We first consider the increment of the two winding numbers when the
radius increases. When $r$ is increased to $r^{\prime}$ ($r^{\prime}>r$),
and the increment $r^{\prime}-r$ is small enough \footnote{Here, ``small enough'' means that the special winding loop picked
in Fig. \ref{fig:strip-winding}(c) can be continuously transported
from $X(E,r)$ to $X(E,r^{\prime})$ while $r$ increases to $r^{\prime}$.}, the increment of $W_{\text{strip}}$ reads, 
\begin{align}
W_{\text{strip}}\left(E,r^{\prime}\right)-W_{\text{strip}}\left(E,r\right) & =\frac{1}{2\pi}\sum_{j=1}^{k}\left(-1\right)^{\tau_{j}}\left(\phi_{j}^{\prime}-\phi_{j}\right),\label{eq:increment-W-strip}
\end{align}
where $\phi_{j}^{\prime}$ and $\phi_{j}$ are the relative phases
of the $j$-th PMGBZ pairs on $X(E,r)$ and $X(E,r^{\prime})$, respectively.
Compared with Eq. (\ref{eq:increment-q1D-winding-number}), the only
difference between the increments of $W$ and $W_{\text{strip}}$
is the sign of $\phi_{j}^{\prime}-\phi_{j}$. In Sec. \ref{app:topo-charge-change-phi}
of the SM \cite{supp-2SK}, we prove that the sign of $\phi_{j}^{\prime}-\phi_{j}$
is related to the topological charge of the corresponding PMGBZ pair.
As illustrated in Fig. \ref{fig:strip-winding}(d), when the arc of
$\phi_{j}$ starts from a negative charge, the relative argument decreases
when the radius increases, shown as the green arrows marked by $\phi_{1}$
and $\phi_{1}^{\prime}$. Otherwise, the relative phase increases,
shown as $\phi_{2}$ and $\phi_{2}^{\prime}$ in Fig. \ref{fig:strip-winding}(d).
Therefore, the increment of $W_{\text{strip}}$ is related to the
increment of $W(E,r)$ by, 
\begin{align}
 & W_{\text{strip}}\left(E,r^{\prime}\right)-W_{\text{strip}}\left(E,r\right)\nonumber \\
= & \frac{1}{2\pi}\sum_{j=1}^{k}\left|\phi_{j}^{\prime}-\phi_{j}\right|,\nonumber \\
= & \frac{W(E,r^{\prime})-W(E,r)}{L_{2}}+O\left(\frac{1}{L_{2}}\right).\label{eq:increment-W-strip-final}
\end{align}

Next, we check the relation between $W_{\text{strip}}$ and $W/L_{2}$
at the limit of $r\rightarrow0^{+}$ and $r\rightarrow+\infty$. In
Sec. \ref{app:relation-degrees} of the SM \cite{supp-2SK}, we prove
that the degrees $M$ and $N$ of $\beta_{1}$ in $F\left(E,\beta_{1}\right)$
are related to the degrees of $\beta_{1}$ in $f\left(E,\beta_{1},\beta_{2}\right)$
by, 
\begin{equation}
M=L_{2}M_{1}+O(1),\quad N=L_{2}N_{1}+O(1),
\end{equation}
where $-M_{1}$ and $N_{1}$ are the lowest and highest degrees of
$\beta_{1}$ in $f\left(E,\beta_{1},\beta_{2}\right)$. For $W\left(E,r\right)$,
when $r\rightarrow0^{+}$, the loop $C_{r}$ encloses no $\beta_{1}$-zeros,
so that, 
\begin{equation}
W\left(E,0^{+}\right)=-M=-M_{1}L_{2}+O(1).
\end{equation}
For the same reason, we also have, 
\begin{equation}
W\left(E,+\infty\right)=N=N_{1}L_{2}+O(1).
\end{equation}
For $W_{\text{strip}}\left(E,r\right)$, when $r\rightarrow0^{+}$,
the term with $\beta_{1}^{-M_{1}}$ dominates in $f\left(E,\beta_{1},\beta_{2}\right)$,
so that all the loop winding numbers $w(\theta_{2})$ equals $-M_{1}$
when $r\rightarrow0^{+}$, and consequently, 
\begin{equation}
W_{\text{strip}}\left(E,0^{+}\right)=\frac{W\left(E,0^{+}\right)}{L_{2}}+O\left(\frac{1}{L_{2}}\right).\label{eq:extreme-case-0}
\end{equation}
Similarly, when $r\rightarrow+\infty$, the term with $\beta_{1}^{N_{1}}$
dominates, so that Eq. (\ref{eq:extreme-case-0}) also holds for $r\rightarrow\infty$.
Combining Eq. (\ref{eq:extreme-case-0}) with Eq. (\ref{eq:increment-W-strip-final}),
we arrive at the conclusion that $W_{\text{strip}}(E,r)$ is the limit
of $W(E,r)/L_{2}$ when $L_{2}\rightarrow\infty$. This relation is
also verified numerically in Sec. \ref{app:num-W-strip} of the SM
\cite{supp-2SK}. Therefore, we get the following SGBZ constraint,
\begin{equation}
\begin{cases}
W_{\text{strip}}\left(E,r_{>}\right)>0\\
W_{\text{strip}}\left(E,r_{<}\right)<0
\end{cases},\label{eq:cond-strip-GBZ}
\end{equation}
which is uniquely determined by the 2D ChP $f(E,\beta_{1},\beta_{2})$.
We define SGBZ as the PMGBZ points on $X(E,r)$, where there exist
$r_{<}\in(r-\epsilon,r)$ and $r_{>}\in(r,r+\epsilon)$ that satisfy
Eq. (\ref{eq:cond-strip-GBZ}) for arbitrary $\epsilon>0$.

Furthermore, the discussion above also shows that $W_{\text{strip}}(E,r)$
increases monotonically with $r$, and has opposite signs when $r\rightarrow0^{+}$
and $r\rightarrow+\infty$. Therefore, the curve of $W_{\text{strip}}(E,r)$
as a function of $r$ will pass the line $W_{\text{strip}}=0$ for
arbitrary reference energy $E$. The difference between the reference
energies out of the SGBZ spectrum and within the SGBZ spectrum lies
in the existence of a platform in the curve of $W_{\text{strip}}(E,r)$
at $W_{\text{strip}}=0$. As illustrated in Fig. \ref{fig:2D-GBZ-condition-and-preimage-of-E}(a),
when the curve exhibits a platform at $W_{\text{strip}}=0$ (orange
region) for some reference energy $E_{1}$, for each point $r_{0}$
on the platform (including the endpoints), $W_{\text{strip}}(E_{1},r)=0$
holds in either $(r_{0}-\epsilon,r_{0})$ or $(r_{0},r_{0}+\epsilon)$
when $\epsilon$ is small enough. Therefore, the SGBZ constraint is
not satisfied and $E_{1}$ does not belong to the SGBZ spectrum, sketched
as the inset of Fig. \ref{fig:2D-GBZ-condition-and-preimage-of-E}(a).
In contrast, when the curve of $W_{\text{strip}}$ does not exhibit
a platform for some reference energy $E_{2}$, as illustrated in Fig.
\ref{fig:2D-GBZ-condition-and-preimage-of-E}(b), the SGBZ constraint
is satisfied at the zero point $r_{0}$ of $W_{\text{strip}}(E_{2},r)$,
and the energy $E_{2}$ belongs to the SGBZ spectrum.

The existence of the zero platform is also manifested in the structure
of $\text{PMGBZ}(E)$. Figures \ref{fig:2D-GBZ-condition-and-preimage-of-E}(c)
and \ref{fig:2D-GBZ-condition-and-preimage-of-E}(d) show the PMGBZ
points with eigenenergy $E_{1}$ and $E_{2}$ on the $\beta_{1}$
complex plane, respectively. For $\text{PMGBZ}(E_{1})$, there is
a circle centered at $0$, which does not intersect with $\text{PMGBZ}(E_{1})$.
Supposing the radius of this circle is $r_{0}$, we define this circle
as a ``central circle'' when $W_{\text{strip}}(E_{1},r_{0})=0$.
Due to the topological invariance of the winding numbers, $W_{\text{strip}}(E_{1},r)$
remains zero in a small neighborhood of $r_{0}$, resulting in a platform
at $W_{\text{strip}}=0$ in the curve of $W_{\text{strip}}(E_{1},r)$.
In contrast, for $\text{PMGBZ}(E_{2})$, there are no central circles,
so the curve of $W_{\text{strip}}(E_{2},r)$ does not show platforms
at $W_{\text{strip}}=0$. Therefore, we arrive at the conclusion that
a reference energy $E$ is out of the SGBZ spectrum if and only if
$\text{PMGBZ}(E)$ exhibits a central circle. The central circle in
our theory shares a deep relation with the central hole in the Amoeba
theory \cite{wangAmoebaFormulationNonBloch2024}, which will be discussed
in detail in Sec. \ref{sec:Amoeba}.

Furthermore, in Sec. \ref{app:winding-general-loop} of the SM \cite{supp-2SK},
we prove that $W_{\text{strip}}(E,r)$ can be calculated by a more
flexible form, 
\begin{equation}
W_{\text{strip}}(E,r)=\frac{1}{2\pi}\int_{-\pi}^{\pi}\tilde{w}\left(s\right)ds,\label{eq:W-strip-new}
\end{equation}
where, 
\begin{equation}
\tilde{w}\left(s\right)=\int_{-\pi}^{\pi}\frac{d\theta_{1}}{2\pi i}\frac{\partial\ln\left[f\left(E,re^{i\theta_{1}},\rho\left(\theta_{1}\right)e^{i\left[s+\delta_{2}\left(\theta_{1}\right)\right]}\right)\right]}{\partial\theta_{1}},\label{eq:def-tilde-w}
\end{equation}
is a winding number around the loop $\theta_{2}=s+\delta_{2}(\theta_{1})$.
Here, $\delta_{2}(\theta_{1}),\theta_{1}\in[-\pi,\pi]$ is an arbitrary
function of $\theta_{1}$ satisfying $\delta_{2}(-\pi)=\delta_{2}(\pi)\mod{2\pi}$.
With the new formulation of Eq. (\ref{eq:def-tilde-w}), $W_{\text{strip}}(E,r)$
can be calculated by winding numbers around arbitrary closed loops
in the form of $\theta_{2}=\delta_{2}(\theta_{1})+s$, which is useful
in some theoretical calculations, such as the calculations in Sec.
\ref{subsec:example}.

\begin{figure}
\centering

\includegraphics{Figures/winding-and-GBZ-candidates-20250426}

\caption{Relation between $\text{PMGBZ}(E)$ and SGBZ spectrum. (a, b) Strip
winding numbers for the reference energy (a) out of the SGBZ spectrum
and (b) inside the SGBZ spectrum. The inset illustrates the relation
between the reference energy and the spectrum. (c, d) The distribution
of (c) $\text{PMGBZ}(E_{1})$ and (d) $\text{PMGBZ}(E_{2})$. The
radius of the red and green dashed circles is $r_{0}$, which is marked
in (a) and (b).}
\label{fig:2D-GBZ-condition-and-preimage-of-E} 
\end{figure}


\section{Anatomy of GDSE with SGBZ formulation}

\label{sec:OBC}

In this section, the relation between the GDSE and the SGBZ is investigated.
Using the 2D Hatano-Nelson (HN) model as an example, we analytically
solve the SGBZs and show that SGBZs can be different in different
strips. Numerical calculations show that both the OBC spectrum and
OBC eigenstates depend on the shapes when the system has different
SGBZs, which implies the existence of the GDSE.

\subsection{Example: 2D HN model with complex coupling coefficients}

\label{subsec:example}

In this part, we will demonstrate the SGBZs of different strips in
a 2D HN model with complex coefficients as an example. The Hamiltonian
of the 2D HN model reads, 
\begin{align}
H_{\text{HN}}= & \sum_{r_{x},r_{y}}J_{x1}c_{r_{x}+1,r_{y}}^{\dagger}c_{r_{x},r_{y}}+J_{x2}c_{r_{x}-1,r_{y}}^{\dagger}c_{r_{x},r_{y}}\nonumber \\
 & +J_{y1}c_{r_{x},r_{y}+1}^{\dagger}c_{r_{x},r_{y}}+J_{y2}c_{r_{x},r_{y}-1}^{\dagger}c_{r_{x},r_{y}},\label{eq:HN-model}
\end{align}
where $J_{x1},J_{x2},J_{y1},J_{y2}$ are complex numbers and $c_{r_{x},r_{y}}$
is the annihilator of the site at coordinate $(r_{x},r_{y})$. The
structure and the coupling terms are illustrated in Fig. \ref{fig:2D-HN-and-stripes}(a).
As shown in Fig. \ref{fig:2D-HN-and-stripes}(b), we choose three
different strips: the $x$-strip with the major axis $\mathbf{a}_{x}=(1,0)$
(purple region), the $y$-strip (green region) with the major axis
$\mathbf{a}_{y}=(0,1)$, and the $[11]$-strip (cyan region) with
the major axis $\mathbf{a}_{[11]}=(1,1)$.

For the $x$-directional and $y$-directional SGBZs, the 2D non-Bloch
ChP under the basis $\{\mathbf{a}_{x},\mathbf{a}_{y}\}$ reads, 
\begin{align}
f_{xy}\left(E,\beta_{x},\beta_{y}\right)= & E-J_{x1}\beta_{x}^{-1}-J_{x2}\beta_{x}\nonumber \\
 & -J_{y1}\beta_{y}^{-1}-J_{y2}\beta_{y},
\end{align}
where $\beta_{x}$ and $\beta_{y}$ are the $\beta$-variables conjugate
to $\mathbf{a}_{x}$ and $\mathbf{a}_{y}$, respectively. We first
construct the base manifold $X(E,r)$. By Vieta's formulas, the two
$\beta_{y}$-zeros satisfy, 
\begin{equation}
\beta_{y}^{(1)}\left(E,\beta_{x}\right)\beta_{y}^{(2)}\left(E,\beta_{x}\right)=\frac{J_{y1}}{J_{y2}},\label{eq:Vieta-theorem-xy}
\end{equation}
where the two solutions are ordered by $|\beta_{y}^{(1)}(E,\beta_{x})|\leq|\beta_{y}^{(2)}(E,\beta_{x})|$.
For any radius $|\beta_{x}|=r$, if we pick $\rho(\theta_{x})=\sqrt{|J_{y1}/J_{y2}|}$,
the relation $|\beta_{y}^{(1)}(E,re^{i\theta_{x}})|\leq\rho(\theta_{x})\leq|\beta_{y}^{(2)}(E,re^{i\theta_{x}})|$
is satisfied, and ``$=$'' holds if and only if the two solutions
have the same modulus. Therefore, $\rho(\theta_{x})=\sqrt{|J_{y1}/J_{y2}|}$
is a valid radius function for $X(E,r)$.

To calculate $w(\theta_{y})$, we consider the value of the ChP restricted
on $X(E,r)$, i.e., 
\begin{align}
f_{xy}\left(\theta_{x},\theta_{y}\right)\equiv & f_{xy}\left(E,re^{i\theta_{x}},\rho(\theta_{x})e^{i\theta_{y}}\right),\nonumber \\
= & E-J_{x1}r^{-1}e^{-i\theta_{x}}-J_{x2}re^{i\theta_{x}}\nonumber \\
 & -J_{y1}\sqrt{\left|\frac{J_{y2}}{J_{y1}}\right|}e^{-i\theta_{y}}-J_{y2}\sqrt{\left|\frac{J_{y1}}{J_{y2}}\right|}e^{i\theta_{y}}.\label{eq:ChP-on-X-xy}
\end{align}
To simplify Eq. (\ref{eq:ChP-on-X-xy}), we make the following substitutions,
\begin{align}
J_{x1} & =\gamma_{x}e^{i\delta_{x}}J_{x},\label{eq:two-types-non-Hermitian-Jx1}\\
J_{x2} & =\gamma_{x}^{-1}e^{i\delta_{x}}J_{x}^{*},\label{eq:two-types-non-Hermitian-Jx2}\\
J_{y1} & =\gamma_{y}e^{i\delta_{y}}J_{y},\label{eq:two-types-non-Hermitian-Jy1}\\
J_{y2} & =\gamma_{y}^{-1}e^{i\delta_{y}}J_{y}^{*},\label{eq:two-types-non-Hermitian-Jy2}
\end{align}
where $\gamma_{x},\gamma_{y}\in\mathbb{R}_{+}$ denote the non-reciprocal
coupling strength, $\delta_{x},\delta_{y}\in\left[-\pi,\pi\right]$
denote the non-reciprocal phases, and $J_{x},J_{y}\in\mathbb{C}$
denote the Hermitian part. With the substitutions above, Eq. (\ref{eq:ChP-on-X-xy})
turns into, 
\begin{align}
f_{xy}\left(\theta_{x},\theta_{y}\right) & =e^{i\delta_{y}}\left[u_{xy}\left(\theta_{y}\right)-v_{xy}\left(\theta_{x}\right)\right],\label{eq:ChP-on-X-simplified-xy}
\end{align}
where, 
\begin{align}
u_{xy}\left(\theta_{y}\right) & =Ee^{-i\delta_{y}}-2\text{Re}\left(J_{y}^{*}e^{i\theta_{y}}\right),\label{eq:uxy}\\
v_{xy}\left(\theta_{x}\right) & =e^{i\Delta_{xy}}\left(\gamma_{x}r^{-1}J_{x}e^{-i\theta_{x}}+\gamma_{x}^{-1}rJ_{x}^{*}e^{i\theta_{x}}\right),\label{eq:vxy}
\end{align}
and $\Delta_{xy}\equiv\delta_{x}-\delta_{y}$. By definition, $w(\theta_{y})$
is the winding number of $f_{xy}$ when $\theta_{y}$ keeps constant
and $\theta_{x}$ winds around $2\pi$, so that it equals the winding
number of $v_{xy}(\theta_{x})$ around $u_{xy}(\theta_{y})$ on the
complex plane. As illustrated in Fig. \ref{fig:2D-HN-and-stripes}(c),
the curve of $u_{xy}$ is a horizontal line segment. The curve of
$v_{xy}$ is an ellipse with its major axis parallel to $\exp(i\Delta_{xy})$.
The length of its major semi-axis is $\gamma_{x}r^{-1}+\gamma_{x}^{-1}r$,
and the length of its minor semi-axis is $|\gamma_{x}^{-1}r-\gamma_{x}r^{-1}|$.
Furthermore, when $r<\gamma_{x}$, $v_{xy}(\theta_{x})$ rotates clockwise
when $\theta_{x}$ increases, and vice versa. Therefore, when $u_{xy}(\theta_{y})$
is enclosed inside the ellipse of $v_{xy}$, shown as the red circle
in Fig. \ref{fig:2D-HN-and-stripes}(c), $w(\theta_{y})=\pm1$, where
the sign of $w(\theta_{y})$ is the same as the sign of $r-\gamma_{x}$.
Otherwise, shown as the green cross in Fig. \ref{fig:2D-HN-and-stripes}(c),
$w(\theta_{y})=0$. As a result, if $u_{xy}$ intersects with $v_{xy}$,
$W_{\text{strip}}(E,r)$ is positive when $r>\gamma_{x}$, and negative
when $r<\gamma_{x}$. According to Eq. (\ref{eq:cond-strip-GBZ}),
$r=\gamma_{x}$ satisfies the SGBZ constraint. The corresponding SGBZ
points are the points on $X(E,\gamma_{x})$ satisfying $u_{xy}(\theta_{y})=v_{xy}(\theta_{x})$.
Thus, we get the SGBZ for the $x$-strip, that is, 
\begin{align}
\beta_{x}\left(\theta_{x},\theta_{y}\right) & =\gamma_{x}e^{i\theta_{x}},\label{eq:SGBZ-xy-x}\\
\beta_{y}\left(\theta_{x},\theta_{y}\right) & =\gamma_{y}e^{i\theta_{y}},\label{eq:SGBZ-xy-y}
\end{align}
where $\theta_{x},\theta_{y}\in[-\pi,\pi]$, and the corresponding
eigenenergy, 
\begin{equation}
E\left(\theta_{x},\theta_{y}\right)=2e^{i\delta_{x}}\text{Re}\left(J_{x}^{*}e^{i\theta_{x}}\right)+2e^{i\delta_{y}}\text{Re}\left(J_{y}^{*}e^{i\theta_{y}}\right).\label{eq:SGBZ-xy-E}
\end{equation}
On the complex plane, the spectrum forms a parallelogram spanned by
$\pm2|J_{x}|e^{i\delta_{x}}$ and $\pm2|J_{y}|e^{i\delta_{y}}$, shown
as Fig. \ref{fig:2D-HN-and-stripes}(d). With the same method, we
can also calculate the SGBZ for the $y$-strip, which is the same
as Eqs. (\ref{eq:SGBZ-xy-x}) and (\ref{eq:SGBZ-xy-y}).

\begin{figure}
\centering

\includegraphics{Figures/HN-example-20250524}

\caption{Schematic diagrams of the 2D HN model, three special strips and the
SGBZ spectra. (a) Illustration of the 2D HN model of Eq. (\ref{eq:HN-model}).
(b) Illustration of the $x$-strip, $y$-strip and $[11]$-strip of
the 2D HN model. (c) Sketches for the calculations of $x$($y$)-SGBZ.
(d) Illustration of the $x$($y$)-SGBZ spectrum. (e) Sketches for
the calculation of $[11]$-SGBZ. (f) Illustration of the $[11]$-SGBZ
spectrum.}
\label{fig:2D-HN-and-stripes} 
\end{figure}

For the $[11]$-directional SGBZ , the 2D non-Bloch ChP under the
basis $\{\mathbf{a}_{[11]},\mathbf{a}_{y}\}$ reads, 
\begin{align}
f_{[11]}\left(E,\beta_{[11]},\beta_{y}\right)= & E-J_{x1}\beta_{[11]}^{-1}\beta_{y}-J_{x2}\beta_{[11]}\beta_{y}^{-1}\nonumber \\
 & -J_{y1}\beta_{y}^{-1}-J_{y2}\beta_{y},
\end{align}
where $\beta_{[11]}$ and $\beta_{y}$ are conjugate to $\mathbf{a}_{[11]}$
and $\mathbf{a}_{y}$, respectively. Similar to the $x$-directional
case, the two $\beta_{y}$-zeros satisfy, 
\begin{equation}
\beta_{y}^{(1)}\left(E,\beta_{[11]}\right)\beta_{y}^{(2)}\left(E,\beta_{[11]}\right)=\frac{J_{x2}\beta_{[11]}+J_{y1}}{J_{x1}\beta_{[11]}^{-1}+J_{y2}}.\label{seq:11-product-beta2-zeros}
\end{equation}
Therefore, the function $\rho(\theta_{[11]})$, defined as, 
\begin{equation}
\rho\left(\theta_{[11]}\right)=\sqrt{\left|\frac{J_{x2}re^{i\theta_{[11]}}+J_{y1}}{J_{x1}r^{-1}e^{-i\theta_{[11]}}+J_{y2}}\right|},\label{eq:rho-theta11-11}
\end{equation}
is a valid radius function for $X(E,r)$. The restricted ChP on $X(E,r)$
reads, 
\begin{align}
f_{[11]}\left(\theta_{[11]},\theta_{y}\right)\equiv & f_{[11]}\left(E,re^{i\theta_{[11]}},\rho(\theta_{[11]})e^{i\theta_{y}}\right),\nonumber \\
= & E-2e^{i\bar{\varphi}(\theta_{[11]})}\sqrt{\left|v_{[11]}\left(\theta_{[11]}\right)\right|}\times\nonumber \\
 & \cos\left(\theta_{y}-\frac{\Delta\varphi(\theta_{[11]})}{2}\right),\label{eq:ChP-on-X-11}
\end{align}
where, 
\begin{align}
v_{[11]}\left(\theta_{[11]}\right)\equiv & J_{x2}J_{y2}re^{i\theta_{[11]}}+J_{x1}J_{y1}r^{-1}e^{-i\theta_{[11]}}\nonumber \\
 & +J_{x1}J_{x2}+J_{y1}J_{y2},
\end{align}
and 
\begin{align}
\varphi_{1}\left(\theta_{[11]}\right) & \equiv\text{Arg}\left(J_{x1}r^{-1}e^{-i\theta_{[11]}}+J_{y2}\right),\\
\varphi_{2}\left(\theta_{[11]}\right) & \equiv\text{Arg}\left(J_{x2}re^{i\theta_{[11]}}+J_{y1}\right),\\
\bar{\varphi}\left(\theta_{[11]}\right) & \equiv\frac{\varphi_{1}\left(\theta_{[11]}\right)+\varphi_{2}\left(\theta_{[11]}\right)}{2},\\
\Delta\varphi\left(\theta_{[11]}\right) & \equiv\varphi_{2}(\theta_{[11]})-\varphi_{1}(\theta_{[11]}).
\end{align}
It is noted that $f_{[11]}$ is continuous in $\theta_{[11]}$ because
$\exp\left[i\bar{\varphi}(\theta_{[11]})\right]$ and $\cos\left(\theta_{y}-\frac{\Delta\varphi(\theta_{[11]})}{2}\right)$
in Eq. (\ref{eq:ChP-on-X-11}) flip their signs simultaneously when
$\theta_{[11]}$ passes the branch cut of $\varphi_{1}$ or $\varphi_{2}$.
To simplify Eq. (\ref{eq:ChP-on-X-11}), we calculate the sum of $w(\theta_{y})$
and $w(\theta_{y}+\pi)$ by the winding number of the following product,
\begin{align}
g\left(\theta_{[11]},\theta_{y}\right)\equiv & f_{[11]}\left(\theta_{[11]},\theta_{y}\right)f_{[11]}\left(\theta_{[11]},\theta_{y}+\pi\right),\nonumber \\
= & 4\cos^{2}\left(\theta_{y}-\frac{\Delta\varphi(\theta_{[11]})}{2}\right)\times\nonumber \\
 & \left[u_{[11]}\left(\theta_{y}-\frac{\Delta\varphi(\theta_{[11]})}{2}\right)-v_{[11]}\left(\theta_{[11]}\right)\right],\label{eq:fun-g}
\end{align}
where $u_{[11]}(s)\equiv E^{2}/4\cos^{2}s$. Therefore, the winding
number of $g$ is related to the winding number of $f$ by, 
\begin{align}
w_{g}\left(\theta_{y}\right)\equiv & \int_{-\pi}^{\pi}\frac{d\theta_{[11]}}{2\pi}\frac{\partial\ln\left[g\left(\theta_{[11]},\theta_{y}\right)\right]}{\partial\theta_{[11]}},\nonumber \\
= & w\left(\theta_{y}\right)+w\left(\theta_{y}+\pi\right),
\end{align}
and the strip winding number equals the integral of $w_{g}(\theta_{y})$
on the interval $[-\pi/2,\pi/2]$.

According to Eq. (\ref{eq:fun-g}), when $\theta_{[11]}$ winds around
$2\pi$, the winding number of $g$ is equal to the winding number
of $v_{[11]}(\theta_{[11]})$ around $u_{[11]}\left(\theta_{y}-\frac{\Delta\varphi(\theta_{[11]})}{2}\right)$.
However, different from the case of the $x$-SGBZ, the argument of
$u_{[11]}$ varies with $\theta_{[11]}$. To fix this problem, we
take the winding loops $\theta_{y}=\Delta\varphi(\theta_{[11]})/2+s$
and $\theta_{y}=\Delta\varphi(\theta_{[11]})/2+s+\pi$ instead of
the loops with constant $\theta_{y}$. Owing to the homology invariance
of the winding numbers (see Sec. \ref{app:winding-general-loop} of
the SM for details \cite{supp-2SK}), $W_{\text{strip}}$ reads, 
\begin{equation}
W_{\text{strip}}\left(E,r\right)=\frac{1}{2\pi}\int_{-\pi/2}^{\pi/2}\tilde{w}_{g}\left(s\right)ds,\label{eq:W-strip-expression-11}
\end{equation}
where, 
\begin{equation}
\tilde{w}_{g}\left(s\right)\equiv\int_{-\pi}^{\pi}\frac{d\theta_{[11]}}{2\pi i}\frac{\partial\ln\left[g\left(\theta_{[11]},\frac{\Delta\varphi(\theta_{[11]})}{2}+s\right)\right]}{\partial\theta_{[11]}},
\end{equation}
is the winding number of $g$ around the loop $\theta_{y}=\Delta\varphi(\theta_{[11]})/2+s$,
which is equal to the winding number of $v_{[11]}(\theta_{[11]})$
around the fixed point $u_{[11]}(s)$. As illustrated in Fig. \ref{fig:2D-HN-and-stripes}(e),
the curve of $v_{[11]}(\theta_{[11]})$ is an ellipse, and the curve
of $u_{[11]}(s)$ is a ray which is collinear with the origin. When
$r<\gamma_{x}\gamma_{y}$, $v_{[11]}(\theta_{[11]})$ rotates clockwise,
and vice versa. When the point $u_{[11]}(s)$ is enclosed by the ellipse,
shown as the red circle in Fig. \ref{fig:2D-HN-and-stripes}(e), $\tilde{w}_{g}(s)=\pm1$,
where the sign is the same as the sign of $r-\gamma_{x}\gamma_{y}$.
Otherwise, when $u_{[11]}(s)$ is out of the ellipse, $\tilde{w}_{g}(s)=0$.
Therefore, according to Eq. (\ref{eq:W-strip-expression-11}), the
SGBZ constraint is satisfied when and only when $r=\gamma_{x}\gamma_{y}$
and $u_{[11]}$ intersects with $v_{[11]}$. The corresponding SGBZ
points are solutions of $f_{[11]}(\theta_{[11]},\theta_{y})=0$ on
$X(E,\gamma_{x}\gamma_{y})$, which reads, 
\begin{align}
\beta_{[11]}\left(\theta_{[11]},\theta_{y}\right) & =\gamma_{x}\gamma_{y}e^{i\theta_{[11]}},\label{eq:SGBZ-11-1}\\
\beta_{y^{\prime}}\left(\theta_{[11]},\theta_{y}\right) & =\gamma_{y}\sqrt{\left|\frac{e^{i\Delta_{xy}}J_{x}^{*}e^{i\theta_{[11]}}+J_{y}}{e^{i\Delta_{xy}}J_{x}e^{-i\theta_{[11]}}+J_{y}^{*}}\right|}e^{i\theta_{y}},\label{eq:SGBZ-11-2}
\end{align}
where $\theta_{[11]},\theta_{y}\in[-\pi,\pi]$. Here, the prime in
the subscript of $\beta_{y^{\prime}}$ is to distinguish it from the
variable $\beta_{y}$ in the $x$-SGBZ or the $y$-SGBZ. The corresponding
eigenenergy is, 
\begin{align}
E\left(\theta_{[11]},\theta_{y}\right)= & 2e^{i\bar{\varphi}(\theta_{[11]})}\sqrt{\left|v_{[11]}\left(\theta_{[11]}\right)\right|}\times\nonumber \\
 & \cos\left(\theta_{y}-\frac{\Delta\varphi(\theta_{[11]})}{2}\right).\label{eq:SGBZ-11-E}
\end{align}
It is noted that the curves $\pm e^{i\bar{\varphi}}\sqrt{|v_{[11]}|}$
are the two square roots of $v_{[11]}$. Therefore, as illustrated
in Fig. \ref{fig:2D-HN-and-stripes}(f), the spectrum of the $[11]$-SGBZ
is the region swept by the line segments connecting the two square
roots of $v_{[11]}(\theta_{[11]})$ for $\theta_{[11]}\in[-\pi,\pi]$.

The above results imply that the SGBZs in different strips are not
necessarily equivalent even in the same non-Hermitian system. For
equivalent SGBZs, under the coordinate transformation, 
\begin{equation}
\begin{pmatrix}\mathbf{a}_{[11]} & \mathbf{a}_{y}\end{pmatrix}=\begin{pmatrix}\mathbf{a}_{x} & \mathbf{a}_{y}\end{pmatrix}\begin{pmatrix}1 & 0\\
1 & 1
\end{pmatrix},
\end{equation}
$\ln\boldsymbol{\beta}\equiv(\ln\beta_{1},\ln\beta_{2})$ should transform
like the momentum $(k_{1},k_{2})$, that is, 
\begin{equation}
\begin{pmatrix}\ln\beta_{[11]} & \ln\beta_{y^{\prime}}\end{pmatrix}=\begin{pmatrix}\ln\beta_{x} & \ln\beta_{y}\end{pmatrix}\begin{pmatrix}1 & 0\\
1 & 1
\end{pmatrix}.\label{eq:coord-trans-HN}
\end{equation}
However, Eq. (\ref{eq:coord-trans-HN}) does not necessarily hold
for the transformation from the $x$($y$)-SGBZ to the $[11]$-SGBZ.
Comparing Eqs. (\ref{eq:SGBZ-xy-x}-\ref{eq:SGBZ-xy-E}) with Eqs.
(\ref{eq:SGBZ-11-1}-\ref{eq:SGBZ-11-E}), the transformation of Eq.
(\ref{eq:coord-trans-HN}) holds if and only if $\sin\Delta_{xy}=0$.
In this case, the real-space Hamiltonian $H_{\text{HN}}$ is similar
to $\exp(i\delta_{x})H_{s}$ by the similarity transformations $|r_{x},r_{y}\rangle\rightarrow\gamma_{x}^{r_{x}}\gamma_{y}^{r_{y}}|r_{x},r_{y}\rangle$,
where 
\begin{equation}
H_{s}=\sum_{r_{x},r_{y}}J_{x}c_{r_{x}+1,r_{y}}^{\dagger}c_{r_{x},r_{y}}\pm J_{y}c_{r_{x},r_{y}+1}^{\dagger}c_{r_{x},r_{y}}+\text{h.c.},\label{eq:H-sym}
\end{equation}
is a Hermitian matrix. The ``$\pm$'' in Eq. (\ref{eq:H-sym}) depends
on whether $\Delta_{xy}=0$ or $\pi$. Therefore, when $\sin\Delta_{xy}=0$,
$H_{\text{HN}}$ shares the same eigenstates with $H_{s}$, which
does not exhibit the GDSE. Otherwise, the system exhibits different
SGBZs in the $x$($y$)-strip and the $[11]$-strip.

The difference between $x$-SGBZ and $[11]$-SGBZ lies in the way
in which the geometry is extended to infinity. For the $x$-SGBZ,
the thermodynamic limit is taken first along $\mathbf{a}_{x}$, then
along $\mathbf{a}_{y}$, while for the $[11]$-SGBZ, the limit is
taken first along $\mathbf{a}_{[11]}$ and then along $\mathbf{a}_{y}$.
The non-equivalence of different SGBZs in the same model implies that
the thermodynamic limit of a 2D or higher-dimensional non-Hermitian
system is dependent on the way the boundary tends to the infinity,
implying the existence of the GDSE.

Moreover, for ChP $f_{[11]}(E,\beta_{[11]},\beta_{y})$, we can also
take $\mathbf{a}_{y}$ as the major axis and $\mathbf{a}_{[11]}$
as the minor axis. In this case, the radius function is, 
\begin{equation}
\rho(\theta_{y})=\sqrt{\left|\frac{J_{x1}}{J_{x2}}\beta_{y}^{2}\right|}=\gamma_{x}r,
\end{equation}
and the value of ChP restricted on $X(E,r)$ reads, 
\begin{align}
f_{[11]}(\theta_{y},\theta_{[11]}) & \equiv f_{[11]}(E,\rho(\theta_{y})e^{i\theta_{[11]}},re^{i\theta_{y}}),\nonumber \\
 & =e^{i\delta_{x}}\left[u_{yx}\left(\theta_{[11]}-\theta_{y}\right)-v_{yx}\left(\theta_{y}\right)\right],
\end{align}
where, 
\begin{align}
u_{yx}\left(s\right) & =Ee^{-i\delta_{x}}-2\text{Re}\left(J_{x}^{*}e^{is}\right),\\
v_{yx}\left(\theta_{y}\right) & =e^{-i\Delta_{xy}}\left(\gamma_{y}r^{-1}J_{y}e^{-i\theta_{y}}+\gamma_{y}^{-1}rJ_{y}^{*}e^{i\theta_{y}}\right).
\end{align}
Compared with Eqs. (\ref{eq:uxy}) and (\ref{eq:vxy}), $u_{yx}$
and $v_{yx}$ are related to $u_{xy}$ and $v_{xy}$ by exchanging
the subscripts $x$ and $y$ for each parameters, respectively. Taking
$\theta_{[11]}=\theta_{y}+s$ as the winding loop, according to Eq.
(\ref{eq:W-strip-new}), $W_{\text{strip}}$ equals the integral of
the winding number, 
\begin{equation}
\tilde{w}(s)=\int_{-\pi}^{\pi}\frac{d\theta_{y}}{2\pi i}\frac{\partial\ln\left[f_{[11]}\left(\theta_{y},\theta_{y}+s\right)\right]}{\partial\theta_{y}},
\end{equation}
which equals the winding number of $v_{yx}(\theta_{y})$ around $u_{yx}(s)$.
By the same reasoning, the SGBZ constraint requires $r=\gamma_{y}$,
and the SGBZ reads $\beta_{[11]}=\gamma_{x}\gamma_{y}e^{i\theta_{[11]}}$
and $\beta_{y^{\prime}}=\gamma_{y}e^{i\theta_{y}}$, which is equivalent
to the $y$-SGBZ of Eqs. (\ref{eq:SGBZ-xy-x}-\ref{eq:SGBZ-xy-E}).
In fact, it is a general conclusion that the selection of the minor
axis does not influence the SGBZ, which will be proved in Sec. \ref{subsec:Equivalence-minor-axes}.

To verify the relation between the SGBZs and the GDSE, in the following
parts, we will show the numerical results of the QMGBZs and OBC eigensystems
of the 2D HN model under different geometries, and compare them with
the analytical solutions of corresponding SGBZs.

\subsection{SGBZs and QMGBZs}

\begin{figure*}
\centering

\includegraphics{Figures/quasi-1D-comparison-20250220}

\caption{Comparison between SGBZs and QMGBZs. (a-c) Theoretical results of
the SGBZ spectra and numerical results of the QMGBZ spectra for (a)
$x$-strip, (b) $y$-strip and (c) $[11]$-strip. (d-f) Complex momenta
in the theoretical results of the SGBZs and the numerical results
of the QMGBZs for (d) $x$-strip, (e) $y$-strip and (f) $[11]$-strip,
where the panels i and ii of each plot show the major and minor $\beta$-variables,
respectively. The parameters of the 2D HN model are $J_{x1}=1+i$,
$J_{x2}=1.5+1.2i$, $J_{y1}=-1+i$ and $J_{y2}=-1.2-0.5i$.}
\label{fig:comparison-quasi-1D} 
\end{figure*}

By definition, the SGBZ is the thermodynamic limit of the QMGBZ. To
verify this, we numerically calculate the QMGBZs of the 2D HN model
in the $x$-strip, the $y$-strip and the $[11]$-strip, then compare
the numerical results with the analytical results given above.

In numerical calculations, the width for each QMGBZ is set to $15$.
The quasi-1D GBZ is solved by first solving the auxiliary GBZ equations
\cite{yangNonHermitianBulkBoundaryCorrespondence2020}, 
\begin{equation}
\begin{cases}
F\left(E,\beta_{1}\right)=0\\
F\left(E,\beta_{1}e^{i\phi}\right)=0
\end{cases},\label{eq:numerical-aGBZ-equation}
\end{equation}
where $E$ and $\beta_{1}$ are the unknowns, and $\phi$ is the relative
phase ranging from $\left[0,\pi\right]$, then checking the 1D GBZ
constraint, i.e., $|\beta_{1}^{(M)}(E)|=|\beta_{1}^{(M+1)}(E)|$.
The numerical solution of Eq. (\ref{eq:numerical-aGBZ-equation})
is given by the Python software package ``phcpy''\cite{verscheldeAlgorithm795PHCpack1999},
which is based on the polynomial homotopy continuation algorithm.
$\phi$ is selected as $\phi=j\pi/N_{\text{\ensuremath{\phi}}},j=0,1,\dots,N_{\phi}$,
where $N_{\phi}=49$.

By solving the equations of Eq. (\ref{eq:numerical-aGBZ-equation}),
the eigenenergy $E$ and the complex momentum in the major axis are
directly calculated. The complex momentum in the minor axis is calculated
by the profile of the non-Bloch waves. For a given pair of $(E,\beta_{1})$,
the non-Bloch wave $\psi_{E,\beta_{1}}$ is calculated by $\left(E-\mathcal{H}\left(\beta_{1}\right)\right)\psi_{E,\beta_{1}}=0$.
We expand the non-Bloch state under the basis of the hybrid space,
i.e., 
\begin{equation}
\psi_{E,\beta_{1}}=\sum_{j}\psi_{E,\beta_{1}}^{(j)}\left|\beta_{1},j\right>,\label{eq:expand-non-Bloch-on-a2}
\end{equation}
where $|\beta_{1},j\rangle\equiv c_{\beta_{1},j}^{\dagger}|0\rangle$
is the state that one particle occupies the site with coordinate $j\mathbf{a}_{2}$
in the supercell. Then, we assume that the non-Bloch state $\psi_{E,\beta_{1}}^{(j)}$
has the form, 
\begin{equation}
\left|\psi_{E,\beta_{1}}^{(j)}\right|^{2}=C\rho^{2j}\cos\left(kx+\varphi\right),\label{eq:fit-supercell}
\end{equation}
where $C$, $\rho$, $k$ and $\varphi$ are parameters for fitting.
By definition, the fit parameter $\rho$ in Eq. (\ref{eq:fit-supercell})
equals $|\beta_{2}|$.

For generality, the coupling terms are four arbitrary complex numbers
selected as $J_{x1}=1+i$, $J_{x2}=1.5+1.2i$, $J_{y1}=-1+i$, and
$J_{y2}=-1.2-0.5i$. Figure \ref{fig:comparison-quasi-1D} illustrates
the comparison between the analytical solutions of the SGBZ and the
numerical solutions of the QMGBZ in the $x$-strip, $y$-strip and
$[11]$-strip. The spectra of the three strips are shown in Fig. \ref{fig:comparison-quasi-1D}(a-c).
In all three strips, the spectra of the QMGBZs (blue dots) fit well
with the SGBZ spectra (orange patches). For the $x$-strip and the
$y$-strip, as shown in Figs. \ref{fig:comparison-quasi-1D}(a) and
\ref{fig:comparison-quasi-1D}(b), the quasi-1D spectra are parallel
line segments to one pair of sides of the SGBZ spectrum. The difference
between the quasi-1D spectra of the $x$-strip and the $y$-strip
lies in the direction of the parallel line segments. For $[11]$-strip,
as shown in Fig. \ref{fig:comparison-quasi-1D}(c), the quasi-1D spectrum
is formed by curves with constant $\theta_{y}-\Delta\varphi(\theta_{[11]})/2$
in Eq. (\ref{eq:SGBZ-11-E}).

Figure \ref{fig:comparison-quasi-1D}(d-f) illustrates the complex
momenta of the QMGBZs and the SGBZs in (d) $x$-strip, (e) $y$-strip
and (f) $[11]$-strip, where the panels i show the moduli of the major
components, and the panels ii show the moduli of the minor components.
For both the $x$-strip and the $y$-strip, the moduli of both $\beta_{x}$
and $\beta_{y}$ are constant, and the constant values are $\left|\beta_{x}\right|=\sqrt{\left|J_{x1}/J_{x2}\right|}\approx0.8580$
and $\left|\beta_{y}\right|=\sqrt{\left|J_{y1}/J_{y2}\right|}\approx1.043$.
For the $[11]$-strip, the modulus of $\beta_{[11]}$ remains constant,
which is $\left|\beta_{[11]}\right|=\sqrt{\left|J_{x1}J_{y1}/J_{x2}J_{y2}\right|}\approx0.8949$,
but the modulus of $\beta_{y}$ depends on $\theta_{[11]}$. For all
three strips, the numerical results fit well with the theoretical
analysis.

In conclusion, in a 2D non-Hermitian lattice, both the spectrum and
the eigenstates of the QMGBZ are consistent with the SGBZ when the
width of the strip is large enough. This conclusion substantiates
our idea that the QMGBZ converges to the SGBZ defined in Sec. \ref{sec:strip-GBZ}
when the width is large enough.

\begin{figure*}
\includegraphics{Figures/OBC-comparison-20241004}

\caption{Comparison between SGBZ spectra and OBC spectra. (a) Spectra of the
parallelogram region spanned by $x$ and $y$ directions, with $L_{x}:L_{y}$
to be (i) $8:1$, (ii) $1:1$ and (iii) $1:8$. The dashed red line
denotes the boundary of the $x$-SGBZ or $y$-SGBZ spectrum. (b) Spectra
of the parallelogram region spanned by $[11]$ and $y$ directions,
with $L_{[11]}:L_{y}$ to be (i) $8:1$, (ii)$1:1$ and (iii) $1:8$.
The dashed red line denotes the boundary of the $[11]$-SGBZ or $y$-SGBZ
spectrum. The shapes of both parallelograms are shown in the inset
of panels i, and the total number of sites is 12800 for each sample.
(c) Proportion of the $[11]$-SGBZ spectrum in the OBC spectrum for
different size and aspect ratios. The number of $[11]$-SGBZ eigenenergies
is counted by the number of eigenenergies enclosed in the theoretical
solution of $[11]$-SGBZ spectrum, and the sign \textquotedblleft$\approx$\textquotedblright{}
in the legend means that the side lengths are integers with the given
aspect ratio and nearest total number, if the exact total number cannot
be reached with the given aspect ratio. The coupling terms of the
2D HN model are the same as Fig. \ref{fig:comparison-quasi-1D}.}
\label{fig:OBC-comparison} 
\end{figure*}


\subsection{SGBZs and OBC eigensystems in parallelogram regions}

According to the discussions above, a 2D non-Hermitian system admits
different SGBZs in different strips. In this part, we show that the
incompatibility of different SGBZs will result in the GDSE of the
OBC spectra and eigenstates.

To illustrate the effect of the different SGBZs in the same model,
we consider the OBC eigensystem in a parallelogram region. Because
a parallelogram is spanned by two axes, by specifying different major
axes, each parallelogram corresponds to two SGBZs. The two SGBZs can
be compatible or incompatible. For example, in the 2D HN model, the
SGBZs of the parallelogram region spanned by the $x$-axis and the
$y$-axis are compatible, but the SGBZs of the parallelogram spanned
by the $[11]$-axis and the $y$-axis are incompatible. 

Figure \ref{fig:OBC-comparison}(a) shows the spectrum of the parallelogram
region spanned by the $x$-axis and the $y$-axis. For numerical calculation,
the total number of sites is fixed to $N_{\text{tot}}=L_{x}\times L_{y}\approx12800$
\footnote{Here, the sign ``$\approx$'' means that the side lengths are selected
as the one with the given aspect and nearest total number to the given
total number if the aspect and the total number are not satisfied
simultaneously. For example, for the case of $N_{\text{tot}}\approx12800$
and $L_{x}:L_{y}=1:1$, the size is $L_{x}=L_{y}=113$, which is the
closest number to $\sqrt{12800}$.}, and the aspect $L_{x}:L_{y}$ is set to $8:1$ (panel i), $1:1$
(panel ii), and $1:8$ (panel iii), where $L_{x}$ and $L_{y}$ denote
the numbers of sites in the $x$-direction and the $y$-direction,
respectively. According to Fig. \ref{fig:OBC-comparison}(a), the
numerical results of the OBC spectra (blue dots) for different spectra
are the same, and consistent with the analytical solution of the $x$-SGBZ
and the $y$-SGBZ (red dashed lines).

For the incompatible case, the OBC spectra of the parallelogram region
spanned by the $[11]$-axis and the $y$-axis are shown in Fig. \ref{fig:OBC-comparison}(b).
The total number of sites is still $N_{\text{tot}}=L_{[11]}\times L_{y}\approx12800$,
and the aspect $L_{[11]}:L_{y}$ is set to $8:1$ (panel i), $1:1$
(panel ii) and $1:8$ (panel iii), where $L_{[11]}$ denotes the number
of sites along the $[11]$-axis. In this case, the OBC spectrum varies
with the aspect ratio. When $L_{[11]}\gg L_{y}$, shown as panel i
of Fig. \ref{fig:OBC-comparison}(b), the OBC spectrum tends to the
spectrum of the $[11]$-SGBZ, and when $L_{y}\gg L_{[11]}$, shown
as panel iii of Fig. \ref{fig:OBC-comparison}(b), the OBC spectrum
tends to the $y$-SGBZ. For the intermediate case, such as $L_{[11]}:L_{y}=1:1$
(panel ii), the OBC spectrum deviates from the spectra of both SGBZs.

Figure \ref{fig:OBC-comparison}(c) shows the proportion of the $[11]$-SGBZ
eigenenergies in the OBC spectrum for different site numbers and aspect
ratios. In numerical calculations, the total number of sites $N_{\text{tot}}=L_{[11]}\times L_{y}$
is set to $3200$, $7200$ and $12800$ ($2:3:4$ in side lengths).
The number of $[11]$-SGBZ eigenenergies ($N_{[11]}$) is calculated
by counting the numerical results of the OBC eigenenergies contained
in the spectrum of the $[11]$-SGBZ, as illustrated by the inset of
Fig. \ref{fig:OBC-comparison}(c). The numerical results show that
the proportion of $[11]$-SGBZ eigenenergies increases rapidly with
the aspect ratio, indicating that the SGBZ describes the limit case
of the OBC eigensystem for aspect ratios far from $1$.

As discussed in the previous part, the SGBZ is the limit of the QMGBZ,
which is a 1D GBZ along the major axis. Hence, the consistency between
the SGBZ with the OBC spectrum is ensured by the 1D non-Bloch band
theory. However, the numerical results tell that the OBC spectrum
deviates from the SGBZ spectrum in incompatible regions when the aspect
ratio is finite, indicating that the growth of the width competes
with the effect of extending the length in the major axis. 

To study the reason for this competition, we numerically calculate
the OBC eigenstate and the corresponding SGBZ mode. For the incompatible
region spanned by the $[11]$-direction and the $y$-direction, we
consider the OBC eigenstates in the parallelogram region with aspect
ratio $8:1$, and the SGBZ modes of the $[11]$-SGBZ. When an OBC
eigenstate $\psi_{\text{OBC}}$ is selected, the SGBZ mode in the
$[11]$-SGBZ can be calculated from the eigenenergy of $\psi_{\text{OBC}}$
according to Eqs. (\ref{eq:SGBZ-11-1}-\ref{eq:SGBZ-11-E}). Figure
\ref{fig:eigenstates-and-GBZ-modes} illustrates the distribution
of the OBC eigenstate $\psi_{\text{OBC}}$ and the SGBZ mode $\psi_{\text{SGBZ}}$
with eigenenergy $1.19+1.23i$. For better visualization, the real-space
coordinates are transformed into the basis of $\{\mathbf{a}_{[11]},\mathbf{a}_{y}\}$,
that is, $\mathbf{r}=r_{1}\mathbf{a}_{[11]}+r_{2}\mathbf{a}_{y}$,
and the modes are scaled by $\psi\left(r_{1},r_{2}\right)\rightarrow\tilde{\psi}\left(r_{1},r_{2}\right)=\gamma_{x}^{-r_{1}}\gamma_{y}^{-r_{1}-r_{2}}\psi\left(r_{1},r_{2}\right)$
to remove the common exponential factors. The real parts of $\tilde{\psi}_{\text{OBC}}\left(r_{1},r_{2}\right)$
(panel i) and $\tilde{\psi}_{\text{SGBZ}}\left(r_{1},r_{2}\right)$
(panel ii) are shown in Fig. \ref{fig:eigenstates-and-GBZ-modes}(a),
and the zoomed-in views of the upper edge (green box marked by `b')
and the lower edge (brown box marked by `c') are shown in Figs. \ref{fig:eigenstates-and-GBZ-modes}(b)
and \ref{fig:eigenstates-and-GBZ-modes}(c), respectively. According
to the numerical results, $\tilde{\psi}_{\text{SGBZ}}$ fits well
with $\tilde{\psi}_{\text{OBC}}$ in the bulk, but deviates from $\tilde{\psi}_{\text{OBC}}$
near the left and right boundaries.

\begin{figure}
\includegraphics{Figures/eigenstates-decomposition-20241022}

\caption{Numerical results of the scaled eigenstate in the parallelogram region
spanned along $[11]$ and $y$ directions. (a) Full view of the real
part of eigenstate under the coordinates $r_{1}$ and $r_{2}$, scaled
by $\gamma_{x}^{-r_{1}}\gamma_{y}^{-r_{1}-r_{2}}$. The OBC eigenstate
and the SGBZ mode are shown in the panel i and panel ii, respectively.
(b, c) Zoomed-in view of the (b) upper edge and (c) lower edge of
the eigenstate, corresponding to the green and brown boxes marked
by `b' and `c'\textbf{ }in (c), respectively. The coupling terms are
the same as Fig. \ref{fig:comparison-quasi-1D}, and the eigenenergy
of the eigenstate is $1.19+1.23i$.}
\label{fig:eigenstates-and-GBZ-modes} 
\end{figure}

Next, we change the width $L_{[11]}$, and compare the deviations
of the SGBZ mode from the OBC eigenstate. To describe the deviations,
we define the residuals $\mathcal{R}_{u}$ and $\mathcal{R}_{l}$
as, 
\begin{align}
\mathcal{R}_{u}\left(r_{1}\right) & =\sum_{r_{2}=\frac{L_{y}}{2}}^{L_{y}}\frac{\left|\tilde{\psi}_{\text{SGBZ}}\left(r_{1},r_{2}\right)-\tilde{\psi}_{\text{OBC}}\left(r_{1},r_{2}\right)\right|^{2}}{\mathcal{N}_{u}},\label{eq:residual-upper}\\
\mathcal{R}_{l}\left(r_{1}\right) & =\sum_{r_{2}=1}^{\frac{L_{y}}{2}}\frac{\left|\tilde{\psi}_{\text{SGBZ}}\left(r_{1},r_{2}\right)-\tilde{\psi}_{\text{OBC}}\left(r_{1},r_{2}\right)\right|^{2}}{\mathcal{N}_{l}},\label{eq:residual-lower}
\end{align}
where $\mathcal{R}_{u}$ and $\mathcal{R}_{l}$ denote the residuals
at the lower edge and the upper edge, respectively, and the normalization
factors are, 
\begin{align}
\mathcal{N}_{u} & =\sum_{r_{1}=1}^{L_{[11]}}\sum_{r_{2}=\frac{L_{y}}{2}}^{L_{y}}\left|\tilde{\psi}_{\text{SGBZ}}\right|^{2}/L_{[11]},\\
\mathcal{N}_{l} & =\sum_{r_{1}=1}^{L_{[11]}}\sum_{r_{2}=1}^{\frac{L_{y}}{2}}\left|\tilde{\psi}_{\text{SGBZ}}\right|^{2}/L_{[11]}.
\end{align}
In numerical calculations, the eigenstates with $L_{y}=30$ and $50$
are compared to the eigenstate with $L_{y}=40$ (i.e., the eigenstate
shown in Fig. \ref{fig:eigenstates-and-GBZ-modes}). For each parallelogram
region, the eigenstate with the closest eigenenergy to $E=1.19+1.23i$
(i.e., the eigenenergy of the eigenstate shown in Fig. \ref{fig:eigenstates-and-GBZ-modes})
is selected. Figures \ref{fig:competition-of-side-lengths}(a) and
\ref{fig:competition-of-side-lengths}(b) show the numerical results
of $\mathcal{R}_{u}$ and $\mathcal{R}_{l}$, respectively. For $L_{y}=30$
and $L_{y}=40$, the residuals decay from the boundaries to the bulk,
which is consistent with the numerical results shown in Fig. \ref{fig:eigenstates-and-GBZ-modes}.
Compared with the case of $L_{y}=40$, the residual in the case of
$L_{y}=30$ decays more rapidly from the boundaries to the bulk. In
contrast, for $L_{y}=50$, the residual is distributed in the bulk,
indicating that the $[11]$-SGBZ mode cannot describe the OBC eigenstate
in this case.

\begin{figure}
\includegraphics{Figures/competition-of-side-lengths-20241112}

\caption{Deviation of SGBZ modes from OBC eigenstates. (a, b) Numerical results
of the residuals (a) $\mathcal{R}_{u}$ and (b) $\mathcal{R}_{l}$
for different $L_{y}$. (c) Quasi-1D picture of the deviation of the
QMGBZ mode from the OBC eigenstate. (d) Illustration of the effect
of the strip width. The blue dots denote the $\beta_{1}$-solutions
corresponding to the auxiliary non-Bloch waves, and the green dot
represents the QMGBZ mode.}
\label{fig:competition-of-side-lengths} 
\end{figure}

Because the SGBZ is equivalent to the QMGBZ when the width is large
enough, the influence of the width can be understood by the quasi-1D
model along the major axis. For a quasi-1D strip with finite width,
supposing that the quasi-1D ChP is $F(E,\beta_{1})$, the number of
OBC equations at the left (right) boundary is equal to $M$ ($N$),
where $-M$ and $N$ are the lowest and highest degrees of $\beta_{1}$
in $F(E,\beta_{1})$, respectively. For each reference energy $E$,
$F(E,\beta_{1})$ has $M+N$ solutions $\beta_{1}^{(j)}(E),j=1,2,\dots,M+N$,
so there exist $M+N$ non-Bloch waves with eigenenergy $E$. However,
the QMGBZ mode only lies in the subspace spanned by the non-Bloch
waves corresponding to $\beta_{1}^{(M)}(E)$ and $\beta_{1}^{(M+1)}(E)$.
Therefore, as shown in Fig. \ref{fig:competition-of-side-lengths}(c),
when the QMGBZ mode itself cannot satisfy all the $M+N$ boundary
equations simultaneously, the non-Bloch waves corresponding to $\beta_{1}^{(j)}(E),j\neq M,M+1$,
dubbed the ``auxiliary non-Bloch waves'', are required to superpose
with the QMGBZ mode to satisfy the OBC. We suppose that the OBC mode
is expanded as the superposition of the non-Bloch waves, 
\begin{align}
\psi_{\text{OBC}}\left(r_{1},r_{2}\right) & =\sum_{j}C_{j}\psi_{\beta_{1}^{(j)}}\left(r_{1},r_{2}\right),\nonumber \\
 & =\sum_{j}C_{j}\left(\beta_{1}^{(j)}\right)^{r_{1}}\phi_{\beta_{1}^{j}}\left(r_{2}\right),
\end{align}
where $\psi_{\beta_{1}^{(j)}}(r_{1},r_{2})\equiv\langle r_{1},r_{2}|\psi_{\beta_{1}^{(j)}}\rangle$
is the real-space function of the non-Bloch wave and $\phi_{\beta_{1}^{j}}(r_{2})\equiv(\beta_{1}^{(j)})^{-1}\psi_{\beta_{1}^{j}}(1,r_{2})$
is the $\beta_{1}$-independent part. Here, the eigenenergy $E$ is
omitted for brevity. To satisfy the $M$ equations at the left boundary,
the first $M-1$ auxiliary non-Bloch waves should be comparable to
the QMGBZ mode at the left boundary, but the last $N-1$ auxiliary
non-Bloch waves should be negligible at the left boundary, that is,
$|C_{j}/C_{M}|$ is comparable to $1$ for $j<M$, and tends to $0$
for $j>M+1$. Similarly, for right $N$ boundary equations, $|C_{j}/C_{M}|\times|\beta_{1}^{(j)}/\beta_{1}^{(M)}|^{L_{1}}$
is comparable to $1$ for $j>M+1$ and tends to $0$ for $j<M$. Therefore,
for finite $L_{1}$, the crosstalk between the left and right boundaries
can be characterized by the ratios $|\beta_{1}^{(M-1)}/\beta_{1}^{(M)}|^{L_{1}}$
and $|\beta_{1}^{(M+1)}/\beta_{1}^{(M)}|^{-L_{1}}$. When the two
ratios tend to $0$, the left (right) auxiliary non-Bloch waves do
not influence the boundary equations at the right (left) boundary.

In 1D lattices, because the difference between $|\beta_{1}^{(M-1)}|$
($|\beta_{1}^{(M+2)}|$) and $|\beta_{1}^{(M)}|$ is generally a nonzero
finite value, both ratios tend to $0$ when $L_{1}$ is large enough.
However, for a quasi-1D strip of a 2D lattice, criticality arises
when the width of the strip tends to infinity. As illustrated in Fig.
\ref{fig:competition-of-side-lengths}(d), both $|\beta_{1}^{(M-1)}|-|\beta_{1}^{(M)}|$
and $|\beta_{1}^{(M+2)}|-|\beta_{1}^{(M)}|$ tend to $0$ with the
order of $1/L_{2}$ when the width $L_{2}\rightarrow\infty$. Supposing,
\begin{equation}
|\beta_{1}^{(M-1)}|=|\beta_{1}^{(M)}|-\frac{\alpha_{L}}{L_{2}}|\beta_{1}^{(M)}|+O\left(\frac{1}{L_{2}^{2}}\right),
\end{equation}
and, 
\begin{equation}
|\beta_{1}^{(M+2)}|=|\beta_{1}^{(M)}|+\frac{\alpha_{R}}{L_{2}}|\beta_{1}^{(M)}|+O\left(\frac{1}{L_{2}^{2}}\right),
\end{equation}
and assuming that $L_{1}$ tends to the infinity by $L_{1}=KL_{2}$,
the crosstalks read, 
\begin{align}
\lim_{L_{2}\rightarrow\infty}\left|\frac{\beta_{1}^{(M-1)}}{\beta_{1}^{(M)}}\right|^{L_{1}} & =\exp\left(-\alpha_{L}K\right),\label{eq:crosstalk-left}\\
\lim_{L_{2}\rightarrow\infty}\left|\frac{\beta_{1}^{(M+2)}}{\beta_{1}^{(M)}}\right|^{-L_{1}} & =\exp\left(-\alpha_{R}K\right),\label{eq:crosstalk-right}
\end{align}
which are determined by the aspect ratio $K$ rather than the size
of the system.

For compatible SGBZs, the SGBZ mode automatically fits the OBC equations
at the boundaries parallel to the minor axis. However, for incompatible
SGBZs, the auxiliary non-Bloch waves matter. When $L_{1}$ increases,
the boundaries parallel to the minor axis move apart, which decreases
the crosstalk. In contrast, when $L_{2}$ increases, the difference
between $|\beta_{1}^{(M-1)}|$ (or $|\beta_{1}^{(M+1)}|$) and $|\beta_{1}^{(M)}|$
decreases, so that the mismatch at the boundaries spreads further
into the bulk, increasing the crosstalk. Due to the competition of
the two effects, the influences of the boundary terms are not negligible
in the bulk even in the thermodynamic limit.

In conclusion, the GDSE is characterized by the incompatibility of
the SGBZs in a given non-Hermitian system. When the system shows different
SGBZs, a parallelogram region with incompatible SGBZs can be construct
by picking the major axes of two different SGBZs. In this case, the
OBC eigensystem depends on the geometry at least in the parallelogram
region. Otherwise, if all the SGBZs of the system are the same, the
eigensystem of an arbitrary region is consistent with the SGBZ eigensystem. 

However, traversing all possible SGBZs in a given system is time-consuming.
A question naturally arises whether the existence of GDSE can be determined
using only one SGBZ. In the next section, we will derive a sufficient
condition for the GDSE, which only requires the information of one
arbitrary SGBZ. 

\section{Basis transformations and conditions for GDSE}

\label{sec:coord-trans}

By definition, the SGBZ is uniquely determined by the 2D non-Bloch
Hamiltonian $h(\beta_{1},\beta_{2})$, or equivalently the 2D non-Bloch
ChP $f(E,\beta_{1},\beta_{2})\equiv\det[E-h(\beta_{1},\beta_{2})]$.
The information of the strip is encoded in the basis $\{\mathbf{a}_{1},\mathbf{a}_{2}\}$
on which the momenta $\beta_{1}$ and $\beta_{2}$ are defined. Therefore,
the relations of different SGBZs are equivalent to the relations of
the ChPs under basis transformations, which enables us to systematically
compare arbitrary two SGBZs in the same system, and derive the general
conditions for the GDSE.

To formulate the problem, we consider two different strips defined
by the basis $\{\mathbf{a}_{1},\mathbf{a}_{2}\}$ and $\{\tilde{\mathbf{a}}_{1},\tilde{\mathbf{a}}_{2}\}$,
where $\mathbf{a}_{1},\tilde{\mathbf{a}}_{1}$ are the major axes
of the two strips and $\mathbf{a}_{2},\tilde{\mathbf{a}}_{2}$ are
the minor axes. Assuming that the transformation of the two bases
reads, 
\begin{equation}
\begin{pmatrix}\tilde{\mathbf{a}}_{1} & \tilde{\mathbf{a}}_{2}\end{pmatrix}=\begin{pmatrix}\mathbf{a}_{1} & \mathbf{a}_{2}\end{pmatrix}\begin{pmatrix}P_{11} & P_{12}\\
P_{21} & P_{22}
\end{pmatrix},\label{eq:transformation-lattice}
\end{equation}
where $P_{ij}\in\mathbb{Z},i,j=1,2$ are the elements of the transform
matrix, the vector $(\ln\beta_{1},\ln\beta_{2})$ is transformed like
the momentum, that is, 
\begin{equation}
\begin{pmatrix}\ln\tilde{\beta}_{1} & \ln\tilde{\beta}_{2}\end{pmatrix}=\begin{pmatrix}\ln\beta_{1} & \ln\beta_{2}\end{pmatrix}\begin{pmatrix}P_{11} & P_{12}\\
P_{21} & P_{22}
\end{pmatrix}.\label{eq:transformation-beta}
\end{equation}
Correspondingly, the ChPs under the two bases satisfy the relation
of Eq. (\ref{eq:transformation-ChP}), 
\begin{align}
f\left(E,\beta_{1},\beta_{2}\right) & =\tilde{f}\left(E,\tilde{\beta}_{1},\tilde{\beta}_{2}\right)\nonumber \\
 & =\tilde{f}\left(E,\beta_{1}^{P_{11}}\beta_{2}^{P_{21}},\beta_{1}^{P_{12}}\beta_{2}^{P_{22}}\right),\label{eq:transformation-ChP}
\end{align}
where $f$ and $\tilde{f}$ are the ChPs under $\{\mathbf{a}_{1},\mathbf{a}_{2}\}$
and $\{\tilde{\mathbf{a}}_{1},\tilde{\mathbf{a}}_{2}\}$, respectively.

In the following part of this section, we first discuss the transformations
that keep the major axis invariant, as illustrated in Fig. \ref{fig:transformation-fixed-major}(a).
We will show that SGBZs are invariant under transformations of this
type, indicating that the SGBZ is uniquely determined by the major
axis. Then, as illustrated in Fig. \ref{fig:transformation-fixed-major}(b),
the transformations that change the major axes are studied. A sufficient
condition for the GDSE is derived by checking whether the SGBZ is
invariant under these transformations.

\subsection{Invariance of SGBZ under minor axis transformations}

\label{subsec:Equivalence-minor-axes}

When the major axis is fixed, the transformation matrix satisfies
$P_{11}=P_{22}=1$ and $P_{21}=0$, and the ChPs satisfy, 
\begin{equation}
f\left(E,\beta_{1},\beta_{2}\right)=\tilde{f}\left(E,\beta_{1},\beta_{1}^{P_{12}}\beta_{2}\right).\label{eq:ChP-transformation-major-axis}
\end{equation}
Supposing that the $\beta_{2}$-zeros of $f$ are $\beta_{2}^{(j)}(E,\beta_{1}),j=1,2,\dots,M_{2}+N_{2}$,
ordered by $\left|\beta_{2}^{(j)}(E,\beta_{1})\right|\leq\left|\beta_{2}^{(k)}(\beta_{1})\right|,\forall j<k$,
in the transformed strip, 
\begin{equation}
\tilde{\beta}_{2}^{(j)}(E,\beta_{1})=\beta_{1}^{P_{12}}\beta_{2}^{(j)}(E,\beta_{1}),
\end{equation}
are also zeros of $\tilde{f}\left(E,\beta_{1},\tilde{\beta}_{2}\right)$,
and the ordering $|\tilde{\beta}_{2}^{(j)}(E,\beta_{1})|\leq|\tilde{\beta}_{2}^{(k)}(E,\beta_{1})|,\forall j<k$
still holds.

For the base manifold $X(E,r)$, by definition, the radius function
$\rho\left(\theta_{1}\right)$ satisfies $\left|\beta_{2}^{(M_{2})}\left(E,re^{i\theta_{1}}\right)\right|\leq\rho\left(\theta_{1}\right)\leq\left|\beta_{2}^{(M_{2}+1)}\left(E,re^{i\theta_{1}}\right)\right|$.
Then, the image of $X(E,r)$ under the coordinate transformation,
i.e., 
\begin{equation}
\tilde{X}(E,r)=\left\{ \left(re^{i\theta_{1}},\tilde{\beta}_{2}\right)\mid\left|\tilde{\beta}_{2}\right|=\rho\left(\theta_{1}\right)r^{P_{12}}\right\} ,
\end{equation}
is also a valid base manifold in the transformed strip. That is, the
radius function $\tilde{\rho}\left(\theta_{1}\right)\equiv r^{P_{12}}\rho\left(\theta_{1}\right)$
satisfies the relation $|\tilde{\beta}_{2}^{(M_{2})}\left(E,re^{i\theta_{1}}\right)|\leq\tilde{\rho}\left(\theta_{1}\right)\leq|\tilde{\beta}_{2}^{(M_{2}+1)}\left(E,re^{i\theta_{1}}\right)|$.

For the strip winding number in the original and transformed strips,
denoted $W_{\text{strip}}(E,r)$ and $\tilde{W}_{\text{strip}}(E,r)$,
we consider the transformation of the winding loop of $w\left(\theta_{2}\right)$
in $X(E,r)$ into the closed loop on $\tilde{X}(E,r)$. As illustrated
in Fig. \ref{fig:transformation-fixed-major}(c) and \ref{fig:transformation-fixed-major}(d),
for the winding loop $\theta_{2}=\theta_{2,0}$ on $X$ with winding
number $w(\theta_{2,0})$, the transformed winding loop reads, 
\begin{equation}
\tilde{\theta}_{2}\left(\theta_{1}\right)=P_{12}\theta_{1}+\theta_{2,0},\label{eq:transformed-theta2}
\end{equation}
and the ChP satisfies, 
\begin{equation}
f\left(E,re^{i\theta_{1}},\rho\left(\theta_{1}\right)e^{i\theta_{2,0}}\right)=\tilde{f}\left(E,re^{i\theta_{1}},\tilde{\rho}\left(\theta_{1}\right)e^{i\tilde{\theta}_{2}\left(\theta_{1}\right)}\right),\label{eq:relation-ChP}
\end{equation}
for arbitrary $\theta_{1}\in\left[-\pi,\pi\right]$. Therefore, in
the transformed strip, the winding number $\tilde{w}(\theta_{2,0})$
around the loop of Eq. (\ref{eq:transformed-theta2}) equals $w(\theta_{2,0})$.
According to the flexible form of the strip winding number {[}Eqs.
(\ref{eq:W-strip-new}) and (\ref{eq:def-tilde-w}){]}, $\tilde{W}_{\text{strip}}(E,r)$
is equal to the integral of $\tilde{w}(\theta_{2,0})$, and consequently
equal to $W_{\text{strip}}(E,r)$ in the original strip.

Because the strip winding numbers $W(E,r)$ and $\tilde{W}(E,r)$
are equal for arbitrary $E$ and $r$, every SGBZ point in the original
strip are transformed into an SGBZ point in the transformed strip.
Therefore, the SGBZs with the same major axis and different minor
axes are compatible to each other.

\begin{figure}
\centering

\includegraphics{Figures/transformation-20250407}

\caption{Basis transformations and conditions for the GDSE. (a, b) Schematic
diagram of the transformations of (a) minor axes and (b) major axes.
(c, d) Unfolded toric base manifolds of (c) original strip and (d)
transformed strip under the minor axis transformation. (e, f) Pairing
of the PMGBZ points under the major axis transformation. The blue
circles denote the PMGBZs.}
\label{fig:transformation-fixed-major} 
\end{figure}


\subsection{Major axis transformations and conditions for GDSE}

Since the SGBZ is independent of the minor axes, without loss of generality,
we consider the transformations keeping the minor axis invariant,
that is, $P_{11}=P_{22}=1$ and $P_{12}=0$. Under this transformation,
the relation of the ChPs reads, 
\begin{equation}
f\left(E,\beta_{1},\beta_{2}\right)=\tilde{f}\left(E,\beta_{1}\beta_{2}^{P_{21}},\beta_{2}\right).\label{eq:ChP-transformation-minor-axis}
\end{equation}
To derive the sufficient condition for the GDSE, or equivalently the
necessary condition for the absence of the GDSE, we require the image
of an SGBZ point under the transformation is also an SGBZ point of
the transformed strip.

According to the definition of SGBZ, an SGBZ point is necessarily
a PMGBZ point. Therefore, when the GDSE is absent, the SGBZ points
in the original strip must be transformed into PMGBZ points in the
transformed strip. For an SGBZ pair $\{(\beta_{1},\beta_{2}),(\beta_{1},\beta_{2}e^{i\phi})\}$,
the transformed points $(\tilde{\beta}_{1},\tilde{\beta}_{2})=(\beta_{1}\beta_{2}^{P_{21}},\beta_{2})$
and $(\tilde{\beta}_{1}^{\prime},\tilde{\beta}_{2}^{\prime})=(\beta_{1}\beta_{2}^{P_{21}}e^{iP_{21}\phi},\beta_{2}e^{i\phi})$
have different $\tilde{\beta}_{1}$ parts. To ensure the transformed
points to be PMGBZ points, for each point $(\beta_{1},\beta_{2})$
in the original SGBZ, as illustrated in Fig. \ref{fig:transformation-fixed-major}(e)
and \ref{fig:transformation-fixed-major}(f), there must exist another
SGBZ point $(\beta_{1}^{\prime},\beta_{2}^{\prime})$ to pair with
$(\beta_{1},\beta_{2})$, such that the transformed points $\left(\tilde{\beta}_{1},\tilde{\beta}_{2}\right)=\left(\beta_{1}\beta_{2}^{P_{21}},\beta_{2}\right)$
and $\left(\tilde{\beta}_{1}^{\prime},\tilde{\beta}_{2}^{\prime}\right)=\left(\beta_{1}^{\prime}\beta_{2}^{\prime P_{21}},\beta_{2}^{\prime}\right)$
are a pair of PMGBZ points with $\tilde{\beta}_{1}=\tilde{\beta}_{1}^{\prime}$,
i.e., 
\begin{align}
\beta_{1}\beta_{2}^{P_{21}} & =\beta_{1}^{\prime}\beta_{2}^{\prime P_{21}},\label{eq:GISE-necessary-1}\\
\left|\beta_{2}\right| & =\left|\beta_{2}^{\prime}\right|,\label{eq:GISE-necessary-2}
\end{align}
Substituting Eq. (\ref{eq:GISE-necessary-2}) into (\ref{eq:GISE-necessary-1}),
the conditions above are simplified into, 
\begin{align}
\left|\beta_{1}\right| & =\left|\beta_{1}^{\prime}\right|,\label{eq:GISE-nece-sim-1}\\
\left|\beta_{2}\right| & =\left|\beta_{2}^{\prime}\right|,\label{eq:GISE-nece-sim-2}\\
\text{Arg}\left(\frac{\beta_{1}}{\beta_{1}^{\prime}}\right) & =\text{Arg}\left[\left(\frac{\beta_{2}^{\prime}}{\beta_{2}}\right)^{P_{21}}\right].\label{eq:GISE-nece-sim-3}
\end{align}
It is noted that both $(\beta_{1},\beta_{2})$ and $(\beta_{1}^{\prime},\beta_{2}^{\prime})$
are located on the original SGBZ, that is, Eqs. (\ref{eq:GISE-nece-sim-1}-\ref{eq:GISE-nece-sim-3})
do not require the information of the transformed SGBZ. Thus, we reach
a necessary condition for the absence of the GDSE. For an arbitrary
SGBZ point $(\beta_{1},\beta_{2})$, and an arbitrary integer $P_{21}$,
there must exist another SGBZ point $(\beta_{1}^{\prime},\beta_{2}^{\prime})$
with the same eigenenergy $E$ to pair with $(\beta_{1},\beta_{2})$,
satisfying Eqs. (\ref{eq:GISE-nece-sim-1}-\ref{eq:GISE-nece-sim-3}).
When $P_{21}$ runs over all integers, the point $(\beta_{1}^{\prime},\beta_{2}^{\prime})$
also changes with $P_{21}$. As a result, an infinite number of SGBZ
points $(\beta_{1}^{\prime},\beta_{2}^{\prime})$ with constant moduli
$|\beta_{1}^{\prime}|,|\beta_{2}^{\prime}|$ and eigenenergy $E$
are required to pair with an SGBZ point $(\beta_{1},\beta_{2})$.

\begin{figure*}
\centering

\includegraphics{Figures/DDS-and-nBDDS-20250624}

\caption{Illustrations of the DDS and non-Bloch DDS, and their relations with
NHSE or GDSE. (a) The case without DDS, where the BZ spectrum is a
1D curve on the complex plane. For some frequency $\omega_{0}$, the
imaginary part of the eigenenergy {[}denoted as $\text{Im}(E)${]}
on the EFC of $\omega_{0}$ is uniform. (b) The case with DDS, where
the EFC of $\omega_{0}$ exhibits non-uniform $\text{Im}(E)$. For
some complex energy $E_{0}$, the number of degenerate momenta on
the BZ (red dots) is finite. (c) The case without non-Bloch DDS, where
the SGBZ spectrum is a 1D curve, and the EFC of some frequency $\omega_{0}$
exhibits uniform $\text{Im}(E)$, $|\beta_{1}|$ and $|\beta_{2}|$
simultaneously. (d) Two cases of the non-Bloch DDS. For the first
case (panel i), the EFC of $\omega_{0}$ exhibits non-uniform $\text{Im}(E)$,
while for the second case (panel ii), the EFC shows uniform $\text{Im}(E)$
but non-uniform $|\beta_{1}|$ or $|\beta_{2}|$. (e) Example for
the non-Bloch DDS with uniform $\text{Im}(E)$. When $|J_{x}|=|J_{y}|$
but $\sin\Delta_{xy}\protect\neq0$, the spectrum of the $[11]$-SGBZ
has zero area (left panel), so that $\text{Im}(E)$ is uniform for
some frequency $\omega_{0}$. However, the modulus $|\beta_{y}|$
is non-uniform on the EFC (right panel), indicating that the non-Bloch
DDS occurs. The parameters for the numerical results in (e) are $\gamma_{x}=1.2$,
$\gamma_{y}=0.5$, $\delta_{x}=\pi/3$, $\delta_{y}=\pi/6$, $J_{x}=\sqrt{2}$
and $J_{y}=1+i$. The eigenfrequency for the EFC is $\omega_{0}=2$.}
\label{fig:DDS} 
\end{figure*}

The physical picture of the necessary condition derived above can
be interpreted as a non-Bloch generalization of the DDS theory \cite{zhangDynamicalDegeneracySplitting2023}.
To understand the concept of the DDS, we first consider the structure
of the degenerate Bloch waves in a 2D Hermitian lattice. Because the
eigenenergies $E(k_{1},k_{2}),(k_{1},k_{2})\in\text{BZ}$ are real,
the degenerate Bloch waves for some eigenenergy $\omega_{0}$ form
a 1D continuum set on the 2D BZ. Next, when non-Hermitian terms are
added to a Hermitian Hamiltonian, the eigenenergies of the BZ are
extended to the complex plane. To find the Bloch waves oscillating
with the eigenfrequency $\omega_{0}$, the equal-frequency contour
(EFC) defined by $\text{Re}[E(k_{1},k_{2})]=\omega_{0}$ is considered.
Different from the Hermitian case, because no constraint is imposed
on the imaginary part of $E$, the Bloch waves on the EFC are not
necessarily degenerate in the sense of complex eigenenergies. 

The degeneracy of the Bloch waves on the EFC can be determined by
the dimensionality of the BZ spectrum around $\text{Re}(E)=\omega_{0}$.
As illustrated in Fig. \ref{fig:DDS}(a), when the line $\text{Re}(E)=\omega_{0}$
(orange dotted line) intersects the BZ spectrum (blue solid curve)
at a point, all the Bloch waves on the EFC share the same imaginary
part of the eigenenergy, determined by the intersection between the
BZ spectrum and the line $\text{Re}(E)=\omega_{0}$. In contrast,
when the intersection is a line segment, shown as Fig. \ref{fig:DDS}(b),
the Bloch waves on the EFC exhibit different imaginary part of eigenenergies.
Therefore, for some fixed complex energy $E_{0}$ (red dot in the
BZ spectrum), the degenerate Bloch waves with eigenenergy $E_{0}$
only form a finite set (red dots on the EFC). The degeneracy splitting
caused by the non-uniform $\text{Im}(E)$ on the EFC is defined as
the DDS. Compared to the degeneracy splitting in Hermitian systems,
these non-degenerate Bloch waves on the EFC oscillates with the same
frequency, but exhibit different decay or growth rates. Therefore,
the superposition of these non-degenerate Bloch waves varies with
time in spite of the coherence of these Bloch waves. According to
Ref. \cite{zhangUniversalNonHermitianSkin2022} and Ref. \cite{zhangDynamicalDegeneracySplitting2023},
the existence of the NHSE is equivalent to the existence of DDS in
non-Hermitian lattices.

With the concept of the DDS, we go back to the condition we derived.
When the system does not exhibit the GDSE, for each SGBZ point $(\beta_{1},\beta_{2})$,
an infinite number of SGBZ points $(\beta_{1}^{\prime},\beta_{2}^{\prime})$
with the same eigenenergy and moduli are required to pair with $(\beta_{1},\beta_{2})$.
Therefore, as illustrated in Fig. \ref{fig:DDS}(c), for each eigenenergy
$E$ with $\text{Re}(E)=\omega_{0}$, there must exist a continuum
of points on the SGBZ with uniform $|\beta_{1}|$, $|\beta_{2}|$
and $\text{Im}(E)$. Because the degenerate SGBZ points form a 1D
continuum on the 2D SGBZ for each eigenenergy, the dimensionality
of the SGBZ spectrum is at most $1$, shown as the left panel of Fig.
\ref{fig:DDS}(c). Therefore, for the EFC on the SGBZ, the absence
of GDSE requires all the non-Bloch waves on the EFC simultaneously
exhibit uniform $\text{Im}(E)$, $|\beta_{1}|$ and $|\beta_{2}|$.
Compared to Fig. \ref{fig:DDS}(a), the absence of the GDSE not only
requires the uniform imaginary parts of the eigenenergies, but also
the uniform distributions of $|\beta_{1}|$ and $|\beta_{2}|$, or
equivalently the imaginary parts of the complex momenta $k_{1}\equiv\ln\beta_{1}$
and $k_{2}\equiv\ln\beta_{2}$. In fact, because the momenta on the
BZ satisfy $\text{Im}k_{1}=\text{Im}k_{2}=0$, Fig. \ref{fig:DDS}(a)
is a special case of Fig. \ref{fig:DDS}(c) when the SGBZ is the same
as the BZ. 

When any one of $\text{Im}(E)$, $|\beta_{1}|$ or $|\beta_{2}|$
is non-uniform on the EFC for some frequency $\omega_{0}$, shown
as Fig. \ref{fig:DDS}(d), Eqs. (\ref{eq:GISE-nece-sim-1}-\ref{eq:GISE-nece-sim-3})
fail, and consequently the GDSE occurs. In this case, we define that
the system exhibits the ``non-Bloch DDS''. Classified by $\text{Im}(E)$
on the EFC, there are two different cases of the non-Bloch DDS. The
first case is shown in the panel i of Fig. \ref{fig:DDS}(d), where
the SGBZ spectrum has non-zero area, and $\text{Im}(E)$ is non-uniform
on the EFC. In this case, the degenerate non-Bloch waves (red dots)
for some complex eigenenergy $E_{0}$ form a finite set. The second
case is shown in the panel ii of Fig. \ref{fig:DDS}(d), where the
area of the SGBZ spectrum is zero, but either $|\beta_{1}|$ or $|\beta_{2}|$
is not constant on the EFC. In this case, although all non-Bloch waves
on the EFC share the same complex eigenenergy, the spatial profile
of these non-Bloch waves cannot match with each other due to the different
spatial decay rates.

The second case of the non-Bloch DDS implies that non-zero spectral
area is sufficient for the GDSE, but zero SGBZ spectral area does
not necessarily rule out the GDSE. For example, the $[11]$-SGBZ of
the 2D HN model has zero spectral area when $|J_{x}|=|J_{y}|$ and
$\sin\Delta_{xy}\neq0$, but according to Sec. \ref{subsec:example},
GDSE occurs when $\sin\Delta_{xy}\neq0$. In Fig. \ref{fig:DDS}(e),
the $[11]$-SGBZ spectrum and the EFC for $\omega_{0}=2$ are numerically
calculated for the 2D HN model with parameters $\gamma_{x}=1.2$,
$\gamma_{y}=0.5$, $\delta_{x}=\pi/3$, $\delta_{y}=\pi/6$, $J_{x}=\sqrt{2}$
and $J_{y}=1+i$. As shown in the left panel, the $[11]$-SGBZ spectrum
consists of two 1D line segments, which has zero spectral area. Therefore,
$\text{Im}(E)$ is uniform on the EFC. By theoretical solutions, we
know that $|\beta_{[11]}|=\gamma_{x}\gamma_{y}$ is also constant
on the EFC, but $|\beta_{y}|$ varies with $\theta_{[11]}$. In the
right panel of Fig. \ref{fig:DDS}(e), the distribution of $|\beta_{y}|$
on the EFC is numerically calculated, which is a non-uniform distribution.

To understand why the continuum degeneracy of the SGBZ points prevents
the GDSE, we recall the discussions in Sec. \ref{sec:OBC} that the
GDSE results from the nonconvergent influence of the auxiliary non-Bloch
waves. For a given eigenenergy, if the dimensionality of the SGBZ
eigenstates is finite, a large number ($\sim L_{2}$) of auxiliary
non-Bloch waves are required to meet the boundary conditions. However,
when the degenerate SGBZ points for some eigenenergy form a continuum,
the number of the SGBZ eigenstates increases with $L_{2}$, so that
the superposition of the SGBZ eigenstates meets most of the boundary
conditions. As a result, the deviation of the OBC eigenstates from
the SGBZ eigenstates is suppressed in this case.

In conclusion, by comparing SGBZs with different major axes, the existence
of the non-Bloch DDS is proved to be a sufficient condition for the
GDSE. We conjecture that the necessity also holds based on the picture
of the boundary-term suppression discussed above. For example, in
the 2D HN model, non-Bloch DDS occurs when and only when $\sin\Delta_{xy}\neq0$,
which is also necessary for a 2D HN model to exhibit the GDSE. 

\section{Relation with Amoeba Formulation}

\label{sec:Amoeba}

\begin{figure*}
\centering

\includegraphics{Figures/Amoeba-20241028}

\caption{Relation between SGBZ spectrum and Amoeba spectrum. (a) Sketch of
the Amoeba spectrum and SGBZ spectra, where $\sigma_{\text{Amoeba}}$
denotes the Amoeba spectrum, and $\sigma_{1},\sigma_{2},\dots$ denote
the SGBZ spectra for different strips. (b) Amoeba spectrum ($\sigma_{\text{Amoeba}}$)
and $[11]$-SGBZ spectrum ($\sigma_{[11]}$) of the 2D HN model, where
$E_{a}=1+i$, $E_{b}=2+2i$ and $E_{c}=3+3i$ are three reference
energies. (c) Amoebas of the three reference energies. (d) PMGBZ distributions
of the three reference energies, where the radius of the green and
red circles is $r=\sqrt{|J_{x1}J_{y1}/J_{x2}J_{y2}|}$. (e) $\beta_{2}$-zeros
with respect to the three reference energies in (a) and the circles
in (d). The red dotted lines denote the valid radius functions $\rho(\theta_{[11]})$.
The parameters of the HN model are the same as Fig. \ref{fig:comparison-quasi-1D}.}
\label{fig:relation-Amoeba} 
\end{figure*}

The Amoeba formulation is a method to construct geometry-independent
GBZs in 2D and higher-dimensional non-Hermitian lattices \cite{wangAmoebaFormulationNonBloch2024}.
However, the relation between the Amoeba formulation and the GDSE
remains unclear. In this section, with the SGBZ description of the
GDSE, we will show that the Amoeba can be viewed as a combination
of all SGBZs in a 2D non-Hermitian lattice, as sketched in Fig. \ref{fig:relation-Amoeba}(a).

Before discussing the spectral relations, we give a brief review of
the Amoeba formulation of 2D and higher-dimensional GBZs. We first
introduce the concept of Amoeba. For a $D$-variate Laurent polynomial
$p(\boldsymbol{\beta})$, the Amoeba of $p$ is defined as, 
\begin{align}
\mathcal{A}_{p}\equiv & \left\{ \ln\left|\boldsymbol{\beta}\right||p\left(\boldsymbol{\beta}\right)=0\right\} ,\label{eq:definition-Amoeba}
\end{align}
where $\ln\left|\boldsymbol{\beta}\right|\equiv\left(\ln\left|\beta_{1}\right|,\ln\left|\beta_{2}\right|,\dots,\ln\left|\beta_{D}\right|\right)$.
Because the non-Bloch ChP $f(E,\boldsymbol{\beta})$ of a $D$-dimensional
non-Hermitian lattice is a $D$-variate Laurent polynomial when a
reference energy $E$ is given as a constant, an Amoeba $\mathcal{A}\left(E\right)\equiv\mathcal{A}_{f\left(E,\cdot\right)}$
is specified by $E$ and $f$. According to the above discussion,
$\mathcal{A}\left(E\right)$ is the set of the logarithms of the moduli
of the possible non-Bloch solutions corresponding to the energy $E$.
To determine the spectrum, an analytic tool called ``Ronkin function''
is defined as, 
\begin{equation}
R_{p}\left(\boldsymbol{\mu}\right)=\int_{\mathbb{T}^{D}}\left(\frac{d\boldsymbol{\theta}}{2\pi}\right)^{D}\ln\left|p\left(e^{\boldsymbol{\mu}+i\boldsymbol{\theta}}\right)\right|,\label{eq:Ronkin-function}
\end{equation}
where $\mathbb{T}^{D}$ is the $d$-dimensional torus when all the
variables $\theta_{1},\theta_{2},\dots,\theta_{D}$ wind around $2\pi$.
By direct calculation, the derivative of $R_{p}\left(\boldsymbol{\mu}\right)$
reads, 
\begin{align}
\nu_{j} & \equiv\frac{\partial R_{p}\left(\boldsymbol{\mu}\right)}{\partial\mu_{j}},\nonumber \\
 & =\text{Re}\int_{T^{d}}\left(\frac{d\theta}{2\pi}\right)^{D}\frac{-i\partial_{\theta_{j}}p\left(e^{\boldsymbol{\mu}+i\boldsymbol{\theta}}\right)}{p\left(e^{\boldsymbol{\mu}+i\boldsymbol{\theta}}\right)},\nonumber \\
 & =\int_{\mathbb{T}^{D-1}}\frac{d\theta_{1}\cdots\widehat{d\theta_{j}}\cdots d\theta_{D}}{\left(2\pi\right)^{D-1}}u_{j}\left(\theta_{1},\dots,\widehat{\theta_{j}},\dots,\theta_{d}\right),\label{eq:derivative-Ronkin}
\end{align}
where the variables with hats mean that they are skipped in the variable
list. The winding number $u_{j}\left(\theta_{1},\dots,\widehat{\theta_{j}},\dots,\theta_{D}\right)$
is defined as, 
\begin{equation}
u_{j}\left(\theta_{1},\dots,\widehat{\theta_{j}},\dots,\theta_{D}\right)\equiv\int_{0}^{2\pi}\frac{d\theta_{j}}{2\pi i}\frac{\partial_{\theta_{j}}p\left(e^{\boldsymbol{\mu}+i\boldsymbol{\theta}}\right)}{p\left(e^{\boldsymbol{\mu}+i\boldsymbol{\theta}}\right)},\label{eq:Amoeba-winding}
\end{equation}
which is the winding number of the ChP when $\theta_{j}$ runs around
$2\pi$ and the other $D-1$ $\theta$-variables remain constant.
For the hole in an Amoeba, i.e., some open set of $\boldsymbol{\mu}\notin\mathcal{A}(E)$
but enclosed by $\mathcal{A}(E)$, owing to the topological invariance
of the winding numbers, $\nu_{j}$ remains constant in the hole. Therefore,
an Amoeba hole can be labeled by the gradient $\left(\nu_{1},\nu_{2},\dots,\nu_{D}\right)$
of the Ronkin function $R_{p}$. When the gradient vanishes, the hole
is called a ``central hole''. The Amoeba theory claims that the
OBC spectrum will tend to Amoeba spectrum $\sigma_{\text{Amoeba}}$
in the thermodynamic limit, where the Amoeba spectrum is defined as
the set of reference energy $E$, satisfying that $\mathcal{A}(E)$
exhibits no central holes.

Under coordinate transformations, it is proved that the Amoeba is
covariant with the coordinates, so that the Amoeba spectrum is unique
for a specific non-Hermitian system. However, when the system exhibits
the GDSE, the thermodynamic limit of the OBC spectrum is not unique,
which contradicts the uniqueness of the Amoeba spectrum. To understand
this contradiction, we consider the relation between the Amoeba spectrum
and the SGBZ. Figure \ref{fig:relation-Amoeba}(b) shows the Amoeba
spectrum $\sigma_{[11]}$ and the $[11]$-SGBZ spectrum $\sigma_{[11]}$
of the 2D HN model, where the parameters are the same as the numerical
calculations in Sec. \ref{sec:OBC}. $E_{a}=1+i$, $E_{b}=2+2i$ and
$E_{c}=3+3i$ are three reference energies. The calculations of the
Amoeba spectrum are available in Sec. \ref{app:calculation-Amoeba}
of the SM \cite{supp-2SK}.

For the Amoeba spectrum, as illustrated in Fig. \ref{fig:relation-Amoeba}(c),
when the reference energies are inside the Amoeba spectrum, such as
$E_{a}$ and $E_{b}$, the Amoeba does not exhibit central holes.
Otherwise, such as $E_{c}$, the Amoeba $\mathcal{A}(E_{c})$ exhibits
a hole containing $(\ln\gamma_{x},\ln\gamma_{y})$. We can prove that
the hole is a central hole. When $\boldsymbol{\mu}=(\ln\gamma_{x},\ln\gamma_{y})$,
or equivalently $|\beta_{x}|=\gamma_{x},|\beta_{y}|=\gamma_{y}$,
the ChP reads, 
\begin{equation}
f_{xy}\left(E,\beta_{x},\beta_{y}\right)=E-2e^{i\delta_{x}}\text{Re}\left(J_{x}^{*}e^{i\theta_{x}}\right)-2e^{i\delta_{y}}\text{Re}\left(J_{y}^{*}e^{i\theta_{y}}\right),
\end{equation}
where $\theta_{j}=\text{Arg}(\beta_{j}),j=x,y$. When $\theta_{x}$
is fixed and $\theta_{y}$ runs over $2\pi$, the trajectory of the
ChP is a line segment parallel to $\exp(i\delta_{y})$, so that the
winding number $u_{y}(\theta_{x})$ is equal to $0$. By the same
reasoning, $u_{x}(\theta_{y})$ also vanishes. Therefore, the hole
in $\mathcal{A}(E_{c})$ is a central hole.

For the SGBZ spectrum, as discussed in Sec. \ref{subsec:Strip-winding-number},
the spectrum is related to the existence of a central circle, i.e.
a circle $C_{r}$ that has no intersections with $\text{PMGBZ}(E)$
and the strip winding number $W_{\text{strip}}(E,r)$ vanishes. Figure
\ref{fig:relation-Amoeba}(d) illustrates $\text{PMGBZ}(E)$ of the
$[11]$-strip on the $\beta_{[11]}$ plane for $E=E_{a}$, $E_{b}$
and $E_{c}$, respectively. When the reference energy lies in the
$[11]$-SGBZ spectrum, such as $E_{a}$, $\text{PMGBZ}(E_{a})$ does
not exhibit central circles. Otherwise, like $\text{PMGBZ}(E_{b})$
and $\text{PMGBZ}(E_{c})$, central circles exist.

To understand why $\mathcal{A}(E_{b})$ does not exhibit a central
hole but $\text{PMGBZ}(E_{b})$ exhibits a central circle, as shown
in Fig. \ref{fig:relation-Amoeba}(e), we consider the zeros $\beta_{y}^{(M_{y})}(E,re^{i\theta_{[11]}})$
and $\beta_{y}^{(M_{y}+1)}(E,re^{i\theta_{[11]}})$ of the ChP $f_{[11]}(E,re^{i\theta_{[11]}},\beta_{y})$,
where $r$ is the radius of the circle shown in Fig. \ref{fig:relation-Amoeba}(d).
For $E=E_{a}$, because $C_{r}$ intersects with $\text{PMGBZ}(E_{a})$,
the curve $|\beta_{y}^{(M_{y})}(E,re^{i\theta_{[11]}})|$ intersects
with $|\beta_{y}^{(M_{y}+1)}(E,re^{i\theta_{[11]}})|$. The intersection
points, marked by the green dots in the left panel of Fig. \ref{fig:relation-Amoeba}(e),
are PMGBZ points located on $X(E,r)$, so that $C_{r}$ is not a central
circle. For $E=E_{c}$, the minimum of $|\beta_{y}^{(M_{y}+1)}|$
is greater than the maximum of $|\beta_{y}^{(M_{y})}|$, so that there
exists a horizontal line that is between the two curves but does not
intersect with them. Suppose that the equation of the line is $|\beta_{y}|=\rho_{0}$,
then the point $(\mu_{[11]},\mu_{y})=(\ln r,\ln\rho_{0})$ is inside
a hole of $\mathcal{A}(E_{c})$, because $f_{[11]}(E_{c},re^{i\theta_{[11]}},\rho_{0}e^{i\theta_{y}})$
does not vanish for arbitrary $\theta_{[11]}$ and $\theta_{y}$.
Furthermore, because the circle $|\beta_{y}|=\rho_{0}$ encloses $M_{2}$
zeros, by Cauchy's argument principle, $u_{y}(\theta_{[11]})$ vanishes,
and consequently $\nu_{y}=0$. By definition, $u_{[11]}(\theta_{y})$
equals the loop winding number, so that $\nu_{[11]}$ equals the strip
winding number $W_{\text{strip}}(E_{c},r)$. Therefore, when $C_{r}$
is a central circle, $\nu_{[11]}=W_{\text{strip}}(E_{c},r)$ vanishes,
and the hole containing $(\ln r,\ln\rho_{0})$ is a central hole.
Conversely, if $(\ln r,\ln\rho_{0})$ is in a central hole, we can
always pick $\rho(\theta_{[11]})=\rho_{0}$ as the radius function,
and the relation $W_{\text{strip}}(E,r)=\nu_{[11]}=0$ holds, indicating
that $C_{r}$ is a central circle.

For $E=E_{b}$, as marked by the green dashed line, the minimum of
$|\beta_{y}^{(M_{y}+1)}|$ is less than the maximum of $|\beta_{y}^{(M_{y})}|$.
Therefore, each horizontal line from $|\beta_{y}|=\min|\beta_{y}^{(M_{y})}|$
to $|\beta_{y}|=\max|\beta_{y}^{(M_{y}+1)}|$ contains at least one
zero of the ChP. Therefore, according to the definition of the Amoeba,
$\mathcal{A}(E_{b})$ does not have holes that intersect the line
$\mu_{[11]}=\ln r$. However, because $|\beta_{y}^{(M_{y}+1)}|$ is
pointwise greater than $|\beta_{y}^{(M_{y})}|$, we can always define
the radius function $\rho(\theta_{[11]})$ such that $X(E_{b},r)$
contains no PMGBZ points, shown as the red dotted curve in the middle
of Fig. \ref{fig:relation-Amoeba}(e). If $W_{\text{strip}}(E_{b},r)=0$,
the circle $C_{r}$ in this case is a central circle, but no central
hole is found in this case.

From this example, we observe that the existence of a central hole
will definitely result in the existence of a central circle, but the
opposite is not necessarily correct. This observation holds for general
2D non-Hermitian systems. If the Amoeba exhibits a central hole that
contains the point $(\mu_{1},\mu_{2})$, then, the line $\ln|\beta_{2}|=\mu_{2}$
must be sandwiched between the curves $|\beta_{2}^{(M_{2})}(E,e^{\mu_{1}+i\theta_{1}})|$
and $|\beta_{2}^{(M_{2}+1)}(E,e^{\mu_{1}+i\theta_{1}})|$ because
$u_{2}(\theta_{1})=0$. Taking $\rho(\theta_{1})=e^{\mu_{2}}$ as
the radius function, the strip winding number satisfies $W_{\text{strip}}(E,e^{\mu_{1}})=\nu_{1}=0$.
Therefore, for each point $(\mu_{1},\mu_{2})$ in the central hole,
the circle $C_{\exp(\mu_{1})}$ is a central circle of $\text{PMGBZ}(E)$.

From the above discussions, we conclude that the SGBZ spectra must
be a subset of the Amoeba spectrum, i.e., 
\begin{equation}
\sigma_{\text{Amoeba}}\supset\cup_{j}\sigma_{j},
\end{equation}
where the spectra $\sigma_{j},j=1,2,\dots$ denote all possible SGBZ
spectra in a given 2D non-Hermitian system. Nevertheless, for the
example 2D HN model, the Amoeba spectrum is the same as the $x$($y$)
SGBZ spectrum, so that the relation, 
\begin{equation}
\sigma_{\text{Amoeba}}=\cup_{j}\sigma_{j},\label{eq:Amoeba-conjecture}
\end{equation}
holds for the 2D HN model with arbitrary complex coupling coefficients.
Therefore, despite the lack of proof, we conjecture that Eq. (\ref{eq:Amoeba-conjecture})
holds for all 2D non-Hermitian lattices. If the conjecture is correct,
the Amoeba spectrum is rigorously equivalent to the combination of
all of the SGBZ spectra.

\section{Conclusion}

In this work, we establish a non-Bloch band theory in 2D non-Hermitian
lattices by proposing the SGBZ formulation. The SGBZ encodes the information
of the geometries, providing a quantitative description of the GDSE
in 2D non-Hermitian lattices. With theoretical and numerical analysis,
we demonstrate that the GDSE is the result of the competition between
incompatible SGBZs. That is, the effects of extending the geometry
along the major axes of two conflicting SGBZs compete with each other.
By comparing the SGBZs in different strips, we find that the dimensionality
of the degenerate SGBZ eigenstates determines whether the GDSE exists
or not. It is proved that the degeneracy splitting from a continuum
set of SGBZ eigenstates to a finite set is a sufficient condition
for the GDSE. Our SGBZ formulation also provides a bridge between
the Amoeba formulation and the GDSE. We show that the Amoeba spectrum
can be understood as a combination of all possible SGBZs. Our SGBZ
formulation also provides a universal guide for the future study of
other important non-Hermitian effects in 2D lattice systems, such
as the non-Hermitian band topology.
\begin{acknowledgments}
This work is supported by the National Key R\textbackslash\&D Program
of China (Grant No. 2023YFA1407600), and the National Natural Science
Foundation of China (NSFC) (Grants No. 12275145, No. 92050110, No.
91736106, No. 11674390, and No. 91836302).
\end{acknowledgments}


\bibliography{2D-skin-effect}

\newpage\onecolumngrid

\setcounter{section}{0} 
\global\long\def\thesection{S\arabic{section}}%
 \setcounter{figure}{0} 
\global\long\def\thefigure{S\arabic{figure}}%
 \setcounter{equation}{0} 
\global\long\def\theequation{S\arabic{equation}}%

% \appendix

\section*{Supplementary Materials}

\section{Density of Quasi-1D Zeros on PMGBZ}

\label{app:DOS-on-1D-GBZ}

In the main text, we have given the relation between the density of
quasi-1D zeros and the relative phases of PMGBZ pairs without proof,
that is, given a reference energy $E$, the number of $\beta_{1}$-zeros
of $F(E,\beta_{1})=0$ located on a segment of PMGBZ is proportional
to the change of relative phases of PMGBZ pairs at both endpoints.
In this section, we will derive the relation with a modified method
of Ref. \cite{yokomizoNonBlochBandTheory2019}. It should be noted
that the reference energy $E$ is viewed as a fixed parameter in this
section.

In general, consider the hybrid Hamiltonian in the following form,
\begin{align}
\mathcal{H}\left(\beta_{1}\right)= & \sum_{\mu,\nu=1}^{m}\sum_{r_{2}=1}^{L_{2}}\mathcal{T}_{0,\mu,\nu}\left(\beta_{1}\right)c_{\beta_{1},r_{2},\mu}^{\dagger}c_{\beta_{1},r_{2},\nu}\nonumber \\
 & +\sum_{t_{2}=1}^{t_{R2}}\sum_{r_{2}=1}^{L_{2}-t_{R2}}\sum_{\mu,\nu=1}^{m}\left[\mathcal{T}_{t_{2},\mu,\nu}\left(\beta_{1}\right)c_{\beta_{1},r_{2}+t_{2},\mu}^{\dagger}c_{\beta_{1},r_{2},\nu}+\mathcal{T}_{-t_{2},\mu,\nu}\left(\beta_{1}\right)c_{\beta_{1},r_{2},\mu}^{\dagger}c_{\beta_{1},r_{2}+t_{2},\nu}\right],\label{seq:quasi-1D-Hamiltonian}
\end{align}
where $\mathcal{T}_{t_{2},\mu,\nu}\left(\beta_{1}\right)$ is the
$\beta_{1}$-dependent coupling coefficient, which is, 
\begin{equation}
\mathcal{T}_{t_{2},\mu,\nu}\left(\beta_{1}\right)=\sum_{t_{1}=-t_{R1}}^{t_{R1}}\mathcal{T}_{t_{1},t_{2},\mu,\nu}\beta_{1}^{-t_{1}}.\label{seq:T-beta1}
\end{equation}
The non-Bloch waves in the minor direction is determined by the 2D
non-Bloch Hamiltonian, 
\begin{equation}
h_{\mu\nu}\left(\beta_{1},\beta_{2}\right)=\sum_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},\mu,\nu}\left(\beta_{1}\right)\beta_{2}^{-t_{2}},\label{seq:2D-non-Bloch-H}
\end{equation}
and the corresponding ChP reads, 
\begin{align}
f\left(E,\beta_{1},\beta_{2}\right) & \equiv\det\left[E-h\left(\beta_{1},\beta_{2}\right)\right],\nonumber \\
 & =\begin{vmatrix}E-\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},1,1}\left(\beta_{1}\right)\beta_{2}^{-t_{2}} & -\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},1,2}\left(\beta_{1}\right)\beta_{2}^{-t_{2}} & \cdots & -\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},1,m}\left(\beta_{1}\right)\beta_{2}^{-t_{2}}\\
-\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},2,1}\left(\beta_{1}\right)\beta_{2}^{-t_{2}} & E-\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},2,2}\left(\beta_{1}\right)\beta_{2}^{-t_{2}} & \cdots & -\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},2,m}\left(\beta_{1}\right)\beta_{2}^{-t_{2}}\\
\vdots & \vdots & \ddots & \vdots\\
-\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},m,1}\left(\beta_{1}\right)\beta_{2}^{-t_{2}} & -\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},m,2}\left(\beta_{1}\right)\beta_{2}^{-t_{2}} & \cdots & E-\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},m,m}\left(\beta_{1}\right)\beta_{2}^{-t_{2}}
\end{vmatrix}.\label{seq:det-2D-ChP}
\end{align}
According to Eq. (\ref{seq:det-2D-ChP}), when $\mathcal{T}_{t_{2},\mu,\nu}\left(\beta_{1}\right)$
are non-zero for any $\mu,\nu=1,2,\dots,m$, the highest and lowest
degrees of $\beta_{2}$ in $f\left(E,\beta_{1},\beta_{2}\right)$
are $-M_{2}=-mt_{R2}$ and $N_{2}=mt_{R2}$, respectively. Therefore,
for given values of $E$ and $\beta_{1}$, there are $M_{2}+N_{2}=2mt_{R2}$
$\beta_{2}$-zeros for the eigenvalue equation $f\left(E,\beta_{1},\beta_{2}\right)=0$.
We sort the zeros by $\left|\beta_{2}^{(1)}\right|\leq\left|\beta_{2}^{(2)}\right|\leq\cdots\leq\left|\beta_{2}^{(2mt_{R2})}\right|$.
The corresponding non-Bloch waves in the reciprocal space can be calculated
by, 
\begin{equation}
h\left(\beta_{1},\beta_{2}^{(j)}\right)\tilde{\phi}^{(j)}=E\tilde{\phi}^{(j)},\label{seq:non-Bloch-wave-reciprocal}
\end{equation}
where $\tilde{\phi}^{(j)}=\left(\tilde{\phi}_{1}^{(j)},\tilde{\phi}_{2}^{(j)},\dots,\tilde{\phi}_{m}^{(j)}\right)\in\mathbb{C}^{m}$,
and the corresponding real-space expression reads, 
\begin{align}
\phi_{\mu}^{(j)}\left(\beta_{1},r_{2}\right) & \equiv\langle\beta_{1},r_{2},\mu|\phi^{(j)}\rangle=\tilde{\phi}_{\mu}^{(j)}\left(\beta_{2}^{(j)}\right)^{r_{2}},
\end{align}
where $\left|\beta_{1},r_{2},\mu\right>\equiv c_{\beta_{1},r_{2},\mu}^{\dagger}\left|0\right>,r_{2}=1,2,\dots,L_{2},\mu=1,2,\dots,m$
is the single-particle basis within a supercell.

Next, we assume that the OBC eigenstate of $\mathcal{H}\left(\beta_{1}\right)$
has the form of Eq. (\ref{seq:formal-OBC-eigenstate}), 
\begin{equation}
\left|\psi\right>=\sum_{j=1}^{2mt_{R2}}C_{j}\left|\phi^{(j)}\right>.\label{seq:formal-OBC-eigenstate}
\end{equation}
Substituting Eqs. (\ref{seq:non-Bloch-wave-reciprocal}-\ref{seq:formal-OBC-eigenstate})
into the eigenvalue equation $\mathcal{H}\left(\beta_{1}\right)\left|\psi\right>=E\left|\psi\right>$,
the equations of $C_{j}$ read, 
\begin{align}
\psi_{\mu}\left(-l\right) & \equiv\sum_{j=1}^{2mt_{R2}}C_{j}\tilde{\phi}_{\mu}^{(j)}\left(\beta_{2}^{(j)}\right)^{-l}=0,\label{seq:eq-C-1}\\
\psi_{\mu}\left(L_{2}+1+l\right) & \equiv\sum_{j=1}^{2mt_{R2}}C_{j}\tilde{\phi}_{\mu}^{(j)}\left(\beta_{2}^{(j)}\right)^{L_{2}+1+l}=0,\label{seq:eq-C-2}
\end{align}
for $l=0,1,\dots,t_{R2}-1$ and $\mu=0,1,\dots,m$. Equations (\ref{seq:eq-C-1})
and (\ref{seq:eq-C-2}) are homogeneous linear equations of $C_{j},j=1,2,\dots,mt_{R2}$,
so that the condition for non-zero solutions of $\left|\psi\right>$
requires the determinant of the coefficients to vanish, i.e., 
\begin{equation}
\begin{vmatrix}\tilde{\phi}^{(1)} & \tilde{\phi}^{(2)} & \cdots & \tilde{\phi}^{(2mt_{R2})}\\
\tilde{\phi}^{(1)}\left(\beta_{2}^{(1)}\right)^{-1} & \tilde{\phi}^{(2)}\left(\beta_{2}^{(2)}\right)^{-1} & \cdots & \tilde{\phi}^{(2mt_{R2})}\left(\beta_{2}^{(2mt_{R2})}\right)^{-1}\\
\vdots & \vdots &  & \vdots\\
\tilde{\phi}^{(1)}\left(\beta_{2}^{(1)}\right)^{-t_{R2}+1} & \tilde{\phi}^{(2)}\left(\beta_{2}^{(2)}\right)^{-t_{R2}+1} & \cdots & \tilde{\phi}^{(2mt_{R2})}\left(\beta_{2}^{(2mt_{R2})}\right)^{-t_{R2}+1}\\
\tilde{\phi}^{(1)}\left(\beta_{2}^{(1)}\right)^{L_{2}+1} & \tilde{\phi}^{(2)}\left(\beta_{2}^{(2)}\right)^{L_{2}+1} & \cdots & \tilde{\phi}^{(2mt_{R2})}\left(\beta_{2}^{(2mt_{R2})}\right)^{L_{2}+1}\\
\tilde{\phi}^{(1)}\left(\beta_{2}^{(1)}\right)^{L_{2}+2} & \tilde{\phi}^{(2)}\left(\beta_{2}^{(2)}\right)^{L_{2}+2} & \cdots & \tilde{\phi}^{(2mt_{R2})}\left(\beta_{2}^{(2mt_{R2})}\right)^{L_{2}+2}\\
\vdots & \vdots &  & \vdots\\
\tilde{\phi}^{(1)}\left(\beta_{2}^{(1)}\right)^{L_{2}+t_{R2}} & \tilde{\phi}^{(2)}\left(\beta_{2}^{(2)}\right)^{L_{2}+t_{R2}} & \cdots & \tilde{\phi}^{(2mt_{R2})}\left(\beta_{2}^{(2mt_{R2})}\right)^{L_{2}+t_{R2}}
\end{vmatrix}=0.\label{seq:det-of-C}
\end{equation}
It is noted that $\tilde{\phi}^{(j)}$ in Eq. (\ref{seq:det-of-C})
is an $m$-vector, so the determinant is a $2mt_{R2}\times2mt_{R2}$
determinant. According to the definition of $\tilde{\phi}^{(j)}$,
$\tilde{\phi}^{(j)}$ is independent of $L_{2}$, so the determinant
in Eq. (\ref{seq:det-of-C}) can be expanded in the general form of
Eq. (\ref{seq:new-eq-C}), 
\begin{equation}
\sum_{1\leq j_{1}<j_{2}<\cdots<j_{M_{2}}\leq2M_{2}}\left(\beta_{2}^{(j_{1})}\beta_{2}^{(j_{2})}\cdots\beta_{2}^{(j_{M_{2}})}\right)^{L_{2}}g_{j_{1}j_{2}\dots j_{M_{2}}}\left(E,\beta_{1},\beta_{2}^{(1)},\dots,\beta_{2}^{(2M_{2})}\right)=0,\label{seq:new-eq-C}
\end{equation}
where the functions $g_{j_{1}j_{2}\dots j_{M_{2}}}\left(E,\beta_{1},\beta_{2}^{(1)},\dots,\beta_{2}^{(2M_{2})}\right)$
are independent of $L_{2}$. When $L_{2}$ tends to infinity, by dividing
both sides of Eq. (\ref{seq:new-eq-C}) by $\left(\beta_{2}^{(M_{2}+1)}\beta_{2}^{(M_{2}+2)}\cdots\beta_{2}^{(2M_{2})}\right)^{L_{2}}$,
we obtain, 
\begin{equation}
g_{M_{2}+1,\dots,2M_{2}}\left(E,\beta_{1},\beta_{2}^{(1)},\dots,\beta_{2}^{(2M_{2})}\right)+\left(\frac{\beta_{2}^{(M_{2})}}{\beta_{2}^{(M_{2}+1)}}\right)^{L_{2}}g_{M_{2},M_{2}+2,\dots,2M_{2}}\left(E,\beta_{1},\beta_{2}^{(1)},\dots,\beta_{2}^{(2M_{2})}\right)+\dots=0,\label{seq:simplified-eq}
\end{equation}
where the omitted terms are less than the order of $\left(\beta_{2}^{(M_{2})}/\beta_{2}^{(M_{2}+1)}\right)^{L_{2}}$.
Note that $\beta_{2}^{(j)},j=1,2,\dots,2M_{2}$ are dependent on $E$
and $\beta_{1}$ by $f\left(E,\beta_{1},\beta_{2}^{(j)}\right)=0$,
when $E$ is fixed, the functions $g_{j_{1}j_{2}\dots j_{M_{2}}}$
are univariate functions of $\beta_{1}$. Therefore, if $\left|\beta_{2}^{(M_{2})}\right|\neq\left|\beta_{2}^{(M_{2}+1)}\right|$,
all the terms in Eq. (\ref{seq:simplified-eq}) vanish except for
$g_{M_{2}+1,\dots,2M_{2}}\left(E,\beta_{1},\beta_{2}^{(1)},\dots,\beta_{2}^{(2M_{2})}\right)$,
so that only a finite number of solutions (that is, independent of
$L_{2}$) of $\beta_{1}$ can be solved from Eq. (\ref{seq:simplified-eq}).
Otherwise, if $\left|\beta_{2}^{(M_{2})}\right|=\left|\beta_{2}^{(M_{2}+1)}\right|$,
the first two terms in Eq. (\ref{seq:simplified-eq}) are preserved,
and the number of solutions increases with the order of $L_{2}$ when
$L_{2}\rightarrow\infty$. This is the reason why the non-Bloch waves
satisfying GBZ constraints account for a major amount of OBC eigenstates.

Given a reference energy $E$, we need to know how these solutions
are distributed on the 1D PMGBZ. Consider the relative argument $\exp\left(i\phi\right)\equiv\beta_{2}^{(M_{2}+1)}/\beta_{2}^{(M_{2})}$.
When $\phi$ is determined, the possible values of $\beta_{1}$ are
also determined. Consequently, all functions $g_{j_{1}j_{2}\dots j_{M_{2}}}$
can be viewed as local functions of $\exp\left(i\phi\right)$. In
the thermodynamic limit, Eq. (\ref{seq:simplified-eq}) reads, 
\begin{equation}
\phi=\frac{i}{L_{2}}\ln\frac{g_{M_{2},M_{2}+2,\dots,2M_{2}}\left(e^{i\phi}\right)}{g_{M_{2}+1,\dots,2M_{2}}\left(e^{i\phi}\right)}+\frac{2n\pi}{L_{2}},n=0,\pm1,\pm2,\dots.\label{seq:phi-eq}
\end{equation}
Next, we define $\tilde{g}\left(e^{i\phi}\right)\equiv i\ln\left[g_{M_{2},M_{2}+2,\dots,2M_{2}}\left(e^{i\phi}\right)/g_{M_{2}+1,\dots,2M_{2}}\left(e^{i\phi}\right)\right]$,
and consider two adjacent solutions, 
\begin{align}
\phi_{1} & =\frac{1}{L_{2}}\tilde{g}\left(e^{i\phi_{1}}\right)+\frac{2n\pi}{L_{2}},\\
\phi_{2} & =\frac{1}{L_{2}}\tilde{g}\left(e^{i\phi_{2}}\right)+\frac{2\left(n+1\right)\pi}{L_{2}},
\end{align}
then the difference of the two solutions reads, 
\begin{align}
\phi_{2}-\phi_{1} & =\frac{2\pi}{L_{2}}+\frac{1}{L_{2}}\frac{d\tilde{g}\left(\phi_{1}\right)}{d\phi}\left(\phi_{2}-\phi_{1}\right)+O\left(\frac{1}{L_{2}^{3}}\right)=\frac{2\pi}{L_{2}}+O\left(\frac{1}{L_{2}^{2}}\right),
\end{align}
where the second equation holds because the order of $\phi_{2}-\phi_{1}$
is $O(1/L_{2})$. Therefore, if the relative phase changes by $\Delta\phi$,
the number of the $\beta_{1}$-solutions reads, 
\begin{align}
N_{\text{q1D}} & =\frac{\left|\Delta\phi\right|}{\frac{2\pi}{L_{2}}+O\left(\frac{1}{L_{2}^{2}}\right)}=\frac{L_{2}}{2\pi}\left|\Delta\phi\right|+O\left(1\right),
\end{align}
which proves the relation mentioned in the main text.

\section{Topological Properties of Loop Winding Numbers and Equivalent Formulation
of $W_{\text{strip}}(E,r)$}

\label{app:winding-general-loop}

\begin{figure}
\centering

\includegraphics{Figures/winding-general-loop-20250527}

\caption{Topological properties and equivalence of winding numbers. (a) Sketch
of the loop winding numbers as a homomorphism between fundamental
groups. (b) Relation between the loop winding numbers $w$ and $\tilde{w}$.
The red and green regions denote the regions where $\tilde{w}(\theta_{2})=w(\theta_{2})+1$
and $\tilde{w}(\theta_{2})=w(\theta_{2})-1$, respectively. (c) Simplification
of the winding number around the loop that winds around $\theta_{2}$-axis.
(d) Homology transformation of an $n$-fold winding loop (blue curve
in i) to $n$ congruent one-fold winding loops (blue solid curve and
purple dashed curve in ii).}
\label{sfig:winding-general-loop} 
\end{figure}

In the main text, we have introduced the strip winding number $W_{\text{strip}}\left(E,r\right)$
defined as the follows, 
\begin{equation}
W_{\text{strip}}\left(E,r\right)=\frac{1}{2\pi}\int_{-\pi}^{\pi}w\left(\theta_{2}\right)d\theta_{2},\label{seq:def-W-strip}
\end{equation}
where the loop winding number $w\left(\theta_{2}\right)$ is defined
as, 
\begin{equation}
w\left(\theta_{2}\right)=\frac{1}{2\pi i}\int_{-\pi}^{\pi}d\theta_{1}\frac{\partial\ln\left[f\left(E,re^{i\theta_{1}},\rho\left(\theta_{1}\right)e^{i\theta_{2}}\right)\right]}{\partial\theta_{1}}.\label{seq:def-W-theta2}
\end{equation}
From a geometric perspective, $w\left(\theta_{2}\right)$ is the winding
number of the 2D ChP $f$ when $\theta_{1}$ runs over $2\pi$ and
$\theta_{2}$ keeps constant. However, in some theoretical derivations,
the exact form of $w\left(\theta_{2}\right)$ is difficult to calculate.
Instead, the winding numbers around some irregular closed loops are
rather simpler. Therefore, a more flexible version of Eq. (\ref{seq:def-W-strip})
is helpful in these cases. In this section, we will show that the
winding number $w\left(\theta_{2}\right)$ in Eq. (\ref{seq:def-W-strip})
can be substituted by $\tilde{w}\left(s\right)$ defined as, 
\begin{equation}
\tilde{w}\left(s\right)=\frac{1}{2\pi i}\int_{-\pi}^{\pi}d\theta_{1}\frac{\partial\ln\left[f\left(E,re^{i\theta_{1}},\rho\left(\theta_{1}\right)e^{i\left[s+\delta_{2}\left(\theta_{1}\right)\right]}\right)\right]}{\partial\theta_{1}},\label{seq:def-tilde-w}
\end{equation}
where $\delta_{2}\left(\theta_{1}\right),\theta_{1}\in\left[-\pi,\pi\right]$
is an arbitrary function of $\theta_{1}$ satisfying $\exp\left[i\delta_{2}\left(-\pi\right)\right]=\exp\left[i\delta_{2}\left(\pi\right)\right]$.

Before discussing the winding number $\tilde{w}\left(s\right)$, we
first give an overview of the topological properties of the winding
number on the toric base space $X(E,r)$. As narrated in the main
text, the toric base space is defined as Eq. (\ref{seq:toric-base-space}),
\begin{equation}
X(E,r)=\left\{ \left(re^{i\theta_{1}},\rho\left(\theta_{1}\right)e^{i\theta_{2}}\right)\in\mathbb{C}^{2}|\theta_{1},\theta_{2}\in\left[-\pi,\pi\right]\right\} ,\label{seq:toric-base-space}
\end{equation}
which is homeomorphic to the 2D torus $\mathbb{T}^{2}$. For a closed
loop $\gamma$ on the space $X$ defined by the parametric equations
$\theta_{1}=\theta_{1}\left(t\right),\theta_{2}=\theta_{2}\left(t\right),t\in\left[0,1\right]$
satisfying $\exp\left[i\theta_{1}\left(0\right)\right]=\exp\left[i\theta_{1}\left(1\right)\right]$
and $\exp\left[i\theta_{2}\left(0\right)\right]=\exp\left[i\theta_{2}\left(1\right)\right]$,
we can define the general loop winding number $w_{\text{loop}}$ defined
as, 
\begin{equation}
w_{\text{loop}}\left(\gamma\right)\equiv\frac{1}{2\pi i}\int_{0}^{1}dt\frac{\partial\ln\left[f\left(E,re^{i\theta_{1}(t)},\rho\left(\theta_{1}(t)\right)e^{i\theta_{2}(t)}\right)\right]}{\partial t}.\label{seq:w-loop}
\end{equation}
The geometric picture of $w_{\text{loop}}$ is illustrated in Fig.
\ref{sfig:winding-general-loop}(a). Because the map $f\left(E,\cdot,\cdot\right):X(E,r)\rightarrow\mathbb{C}$
is a continuous map, closed loops {[}loops $\gamma_{1}$ and $\gamma_{2}$
in Fig. \ref{sfig:winding-general-loop}(a) {]} on $X(E,r)$ are mapped
to closed loops on $\mathbb{C}$ {[}purple and orange closed loops
in the complex plane of Fig. \ref{sfig:winding-general-loop}(a) {]}.
When the image of $\gamma$ under $f\left(E,\cdot,\cdot\right)$ does
not pass $0$, $w_{\text{loop}}\left(\gamma\right)$ is well defined,
and equal to the number of laps that the image of $\gamma$ winds
around $0$. Because the preimage of $0$ under $f\left(E,\cdot,\cdot\right)$
is $\text{PMGBZ}\left(E\right)\cap X(E,r)$, we consider the restriction
of $f\left(E,\cdot,\cdot\right)$ on $X(E,r)\backslash\text{PMGBZ}\left(E\right)$,
that is, $f\left(E,\cdot,\cdot\right):X(E,r)\backslash\text{PMGBZ}\left(E\right)\rightarrow\mathbb{C}\backslash\left\{ 0\right\} $.
When $\text{PMGBZ}\left(E\right)\cap X(E,r)$ is a finite set, $X(E,r)\backslash\text{PMGBZ}\left(E\right)$
is pathwise connected, so that $f\left(E,\cdot,\cdot\right)$ induces
a homomorphism $f_{*}\left(E,\cdot,\cdot\right):\pi_{1}\left(X(E,r)\backslash\text{PMGBZ}\left(E\right)\right)\rightarrow\pi_{1}\left(\mathbb{C}\backslash\left\{ 0\right\} \right)=\pi_{1}\left(\mathbb{T}^{1}\right)=\mathbb{Z}$,
which maps the homotopy class $\left[\gamma\right]$ to $w_{\text{loop}}\left(\gamma\right)$
\cite{hatcherAlgebraicTopology2001}. In other words, the map $w_{\text{loop}}:\left\{ \text{Closed loops on }X(E,r)\backslash\text{PMGBZ}\left(E\right)\right\} \rightarrow\mathbb{Z}$
satisfies $w_{\text{loop}}\left(\gamma_{1}\circ\gamma_{2}\right)=w_{\text{loop}}\left(\gamma_{1}\right)+w_{\text{loop}}\left(\gamma_{2}\right)$,
where $\circ$ is the product of the fundamental group $\pi_{1}\left(X(E,r)\backslash\text{PMGBZ}\left(E\right)\right)$.
Furthermore, because $\mathbb{Z}$ is abelian, $w_{\text{loop}}$
vanishes on the commutators, that is, $w_{\text{loop}}\left(\gamma_{1}\circ\gamma_{2}\circ\gamma_{1}^{-1}\circ\gamma_{2}^{-1}\right)=0,\forall\gamma_{1},\gamma_{2}\in X(E,r)\backslash\text{PMGBZ}(E)$,
so that $w_{\text{loop}}$ is also a homomorphism from the homology
group $H_{1}\left(X(E,r)\backslash\text{mGBZ}\left(E\right)\right)$
\cite{hatcherAlgebraicTopology2001} to the integer group $\mathbb{Z}$.

The topological nature of $w_{\text{loop}}$ allows us to perform
algebraic operations on the winding numbers of the closed loops. As
an example, the difference $\gamma_{2}-\gamma_{1}$ (in the sense
of homology) in Fig. \ref{sfig:winding-general-loop}(a) is homologous
to an infinitesimal loop around the PMGBZ point, shown as the red
arrowed loop in Fig. \ref{sfig:winding-general-loop}(a). As a result,
the difference $w_{\text{loop}}\left(\gamma_{2}\right)-w_{\text{loop}}\left(\gamma_{1}\right)$
is equal to the winding number around the infinitesimal loop, which
is defined as the topological charge of the PMGBZ point.

Based on the discussion above, we will show that $W_{\text{strip}}$
can also be calculated by the following integral, 
\begin{equation}
W_{\text{strip}}\left(E,r\right)=\frac{1}{2\pi}\int_{-\pi}^{\pi}\tilde{w}\left(s\right)ds.\label{seq:W-strip-new}
\end{equation}
As defined in Eq. (\ref{seq:def-tilde-w}), for each constant value
$s$, $\tilde{w}\left(s\right)$ equals the winding number of the
loop defined by $\theta_{2}\left(\theta_{1}\right)=\delta_{2}\left(\theta_{1}\right)+s$
and $\theta_{1}=2\pi t,t\in\left[0,1\right]$. We first consider the
case where $\delta_{2}\left(-\pi\right)=\delta_{2}\left(\pi\right)$,
that is, the loop does not wind around the $\theta_{2}$-axis. Without
loss of generality, we suppose $\delta_{2}(-\pi)=\delta_{2}(\pi)=0$.
As shown in Fig. \ref{sfig:winding-general-loop}(b), for each constant
value $s$, the difference between $\tilde{w}\left(s\right)$ and
$w\left(s\right)$ equals the total topological charges enclosed by
the loop of $\tilde{w}\left(s\right)$ {[}orange solid line in Fig.
\ref{sfig:winding-general-loop}(b){]} and the reverse of the loop
of $w\left(s\right)$ {[}orange dashed line in Fig. \ref{sfig:winding-general-loop}(b){]}.
For a pair of PMGBZ points, as illustrated by the light red and light
green regions in Fig. \ref{sfig:winding-general-loop}(b), the contributions
of the positive charge and the negative charge to $\int_{-\pi}^{\pi}\left[\tilde{w}\left(s\right)-w\left(s\right)\right]ds$
cancel with each other, so that the integral of $\tilde{w}\left(s\right)$
equals the integral of $w\left(s\right)$, which proves Eq. (\ref{seq:W-strip-new}).

Next, we consider the general case where the loop defined by $\theta_{2}=\delta_{2}\left(\theta_{1}\right)$
is allowed to wind around the $\theta_{2}$-axis, shown as the blue
solid line in Fig. \ref{sfig:winding-general-loop}(c). As illustrated
in the figure, the blue loop is homotopic to the purple dashed loop,
and the purple dashed loop is homologous to the sum of the orange
dotted loop, which satisfies the condition $\delta_{2}\left(-\pi\right)=\delta_{2}\left(\pi\right)$,
and the gray loop, around which the winding number of ChP vanishes.
Therefore, for each winding loop (solid blue line), the winding number
around the loop equals the winding number of a special loop satisfying
$\delta_{2}\left(-\pi\right)=\delta_{2}\left(\pi\right)$ (dotted
orange line), which is the case shown in Fig. \ref{sfig:winding-general-loop}(b).

Furthermore, in some cases, the form of the ChP is complicated, but
the product, 
\begin{equation}
g\left(\theta_{1},\theta_{2}\right)=\prod_{m=1}^{n}f\left(E,re^{i\theta_{1}},\rho(\theta_{1})e^{i\left(\theta_{2}+2m\pi/n\right)}\right),
\end{equation}
has a simple form, such as the $[11]$-SGBZ in 2D HN model discussed
in the main text. Then, $W_{\text{strip}}$ can be calculated by,
\begin{equation}
W_{\text{strip}}\left(E,r\right)=\frac{1}{2\pi}\int_{-\pi/n}^{\pi/n}\tilde{w}_{g}\left(s\right)ds,\label{eq:W-strip-g}
\end{equation}
where, 
\begin{equation}
\tilde{w}_{g}\left(s\right)=\int_{-\pi}^{\pi}\frac{d\theta_{1}}{2\pi i}\frac{\partial\ln\left[g\left(\theta_{1},s+\delta_{2}(\theta_{1})/n\right)\right]}{\partial\theta_{1}},\label{seq:w-g}
\end{equation}
is the winding number of $g$ on the loop $s+\delta_{2}(\theta_{1})/n$.
Here, we only require $\delta_{2}(\theta_{1})$ rather than $\delta_{2}(\theta_{1})/n$
to satisfy the periodicity condition $\exp[i\delta_{2}(-\pi)]=\exp[i\delta_{2}(\pi)]$.
Now, we will prove Eq. (\ref{seq:w-g}).

First, $g(\theta_{1},s+\delta_{2}(\theta_{1})/n)$ is a periodic function
of $\theta_{1}$. When $\theta_{1}$ increases by $2\pi$, $\delta_{2}(\theta_{1})/n$
increases by an integer multiple of $2\pi/n$. Supposing $\delta_{2}(\pi)/n=\delta_{2}(-\pi)/n+2\pi m_{0}/n\mod{2\pi}$,
we get, 
\begin{align}
g\left(\pi,s+\delta_{2}(\pi)/n\right) & =\prod_{m=1}^{n}f\left(E,re^{i\pi},\rho(\pi)e^{i\left(s+\delta_{2}(\pi)/n+2m\pi/n\right)}\right),\nonumber \\
 & =\prod_{m=1}^{n}f\left(E,re^{-i\pi},\rho(-\pi)e^{i\left(s+\delta_{2}(-\pi)/n+2(m+m_{0})\pi/n\right)}\right),\nonumber \\
 & =g\left(-\pi,s+\delta_{2}(-\pi)/n\right).
\end{align}
Therefore, the number $\tilde{w}_{g}(s)$ defined by Eq. (\ref{seq:w-g})
is a valid winding number.

Next, we consider the relation between $\tilde{w}_{g}(s)$ and the
winding number of the ChP. When $\delta_{2}(-\pi)/n=\delta_{2}(\pi)/n\mod{2\pi}$,
the curve $\theta_{2}=s+\delta_{2}(\theta_{1})/n$ itself forms a
closed loop on $X(E,r)$. Therefore, the following relation holds,
\begin{align}
\tilde{w}_{g}\left(s\right) & =\sum_{m=1}^{n}\int_{-\pi}^{\pi}\frac{d\theta_{1}}{2\pi i}\frac{\partial\ln\left[f\left(E,re^{i\theta_{1}},\rho(\theta_{1})e^{i\left(s+\delta_{2}(\theta_{1})/n+2m\pi/n\right)}\right)\right]}{\partial\theta_{1}},\nonumber \\
 & =\sum_{m=1}^{n}\tilde{w}\left(s+2m\pi/n\right),\label{seq:relation-winding-numbers}
\end{align}
where $\tilde{w}(s)$ here is defined by the loop $\delta_{2}^{\prime}(\theta_{1})\equiv\delta_{2}(\theta_{1})/n$.
By Eq. (\ref{seq:W-strip-new}), the integral of Eq. (\ref{eq:W-strip-g})
reads, 
\begin{align}
\frac{1}{2\pi}\int_{-\pi/n}^{\pi/n}\tilde{w}_{g}\left(s\right)ds= & \sum_{m=1}^{n}\frac{1}{2\pi}\int_{-\pi/n}^{\pi/n}\tilde{w}\left(s+2m\pi/n\right)ds,\nonumber \\
= & \frac{1}{2\pi}\int_{\pi/n}^{2\pi+\pi/n}\tilde{w}\left(s\right)ds,
\end{align}
which equals $W_{\text{strip}}(E,r)$.

When $\delta_{2}(-\pi)/n=\delta_{2}(\pi)/n+2\pi m_{0}/n\mod{2\pi},m_{0}\neq0$,
the curve $\theta_{2}=s+\delta_{2}(\theta_{1})/n,\theta_{1}\in[-\pi,\pi]$
is not closed, so that the winding number $\tilde{w}(s)$ in Eq. (\ref{seq:relation-winding-numbers})
is ill-defined. However, by homology transformations, $\tilde{w}_{g}(s)$
can be transformed into the sum of the winding numbers around $n$
congruent winding loops with a phase shift of $2\pi/n$ along the
$\theta_{2}$-axis. Taking $n=2,m_{0}=1$ as an example, as shown
in panel i of Fig. \ref{sfig:winding-general-loop}(d), the curve
$\delta_{2}(\theta_{1})/2$ is braided with $\delta_{2}(\theta_{1})/2+\pi$
(black dotted curves), forming a twofold winding loop around the $\theta_{1}$-axis
(solid blue curves). The winding number $\tilde{w}_{g}(s)$ is equal
to the winding number of the ChP around the twofold loop. By virtue
of the homology invariance of the winding number, a winding loop parallel
to the $\theta_{2}$-axis can be added to the two-fold loop without
changing the winding number, shown as the gray horizontal line in
panel i of Fig. \ref{sfig:winding-general-loop} (d). Then, the sum
of the two-fold loop and the horizontal loop is homologous to the
sum of two one-fold loops, shown as the blue solid curve and the purple
dashed curve in panel ii of Fig. \ref{sfig:winding-general-loop}(d).
Thus, we return to the case with $m_{0}=0$. In general, with the
same method, an arbitrary $n$-fold loop can be split into $n$ one-fold
loops by homology, so Eq. (\ref{eq:W-strip-g}) holds for arbitrary
$n$ and $m_{0}$.

\section{Sign of Topological Charge and Change of Relative Argument}

\label{app:topo-charge-change-phi}

\begin{figure}
\centering

\includegraphics{Figures/supp-derivative-zeros-20250120}

\caption{Sketch for the relation between topological charge and relative phase.
(a) 3D view of a PMGBZ pair and the $\beta_{2}$-zeros in the $\theta_{1}$-$\beta_{2}$
space. The gray surface is the toric base manifold $X(E,r)$, $p_{\mathrm{a}}$
and $p_{\mathrm{b}}$ form a PMGBZ pair, and $\beta_{2}^{(\mathrm{a})},\beta_{2}^{(\mathrm{b})}$
are the $\beta_{2}$-zeros of $f(E,re^{i\theta_{1}},\beta_{2}^{(j)})$
in the neighborhoods of $p_{\mathrm{a}}$ and $p_{\mathrm{b}}$, respectively.
(b) Absolute value of $\beta_{2}^{(\mathrm{a})}$ and $\beta_{2}^{(\mathrm{b})}$
in the neighborhood of a PMGBZ pair.}
\label{sfig:topo-charge-rel-arg} 
\end{figure}

In the main text, we show that the increment of strip winding number
is related to the change of relative arguments by, 
\begin{equation}
W_{\text{strip}}\left(E,r^{\prime}\right)-W_{\text{strip}}\left(E,r\right)=\frac{1}{2\pi}\sum_{j=1}^{k}\left(-1\right)^{\tau_{j}}\left(\phi_{j}^{\prime}-\phi_{j}\right),\label{seq:W-strip}
\end{equation}
where $\phi_{j}^{\prime}$ and $\phi_{j}$ are the relative arguments
of the $j$-th PMGBZ pair at radius $r$ and $r^{\prime}$, and $\tau_{j}=0,1$
is related to the topological charge of the PMGBZ pair. If the base
point of the relative phase $\phi_{j}$ takes a positive charge, $\tau_{j}=0$,
and vice versa. In this section, we will show that the sign of $\phi_{j}^{\prime}-\phi_{j}$
is related to the topological charge when $r^{\prime}>r$.

For an PMGBZ pair consisting of two PMGBZ points $p_{\mathrm{a}}=\left(\mathcal{B}_{1},\mathcal{B}_{2}\right)$
and $p_{\mathrm{b}}=\left(\mathcal{B}_{1},\mathcal{B}_{2}e^{i\phi}\right)$,
as schematically illustrated in Fig. \ref{sfig:topo-charge-rel-arg}(a),
there must exist two curves of $\beta_{2}$-zeros of 2D ChP in the
two neighborhoods of $p_{\mathrm{a}}$ and $p_{\mathrm{b}}$, shown
as the curves $\beta_{2}^{(\mathrm{a})}$ and $\beta_{2}^{(\mathrm{b})}$
in Fig. \ref{sfig:topo-charge-rel-arg}(a). It is noted that the ordering
$|\beta_{2}^{(1)}|\leq|\beta_{2}^{(2)}|\leq\cdots\leq|\beta_{2}^{(M_{2}+N_{2})}|$
used in the main text will lead to the discontinuity of $\beta_{2}^{(j)}$
as a function of $\theta_{1}$. For example, the $M_{2}$-th solutions
{[}blue and brown dashed curves in Fig. \ref{sfig:topo-charge-rel-arg}(a){]}
jump into $\left(M_{2}+1\right)$-th solutions {[}blue and brown solid
curves in Fig. \ref{sfig:topo-charge-rel-arg}(a){]} when passing
an PMGBZ point. Therefore, in this section, we use the superscript
$(\mathrm{a})$ or $(\mathrm{b})$ to label the continuous $\beta_{2}$-zero
curve that passes the PMGBZ point $p_{\mathrm{a}}$ or $p_{\mathrm{b}}$,
instead of the indices ordered by absolute values.

The topological charge of a PMGBZ point is determined by the direction
in which the $\beta_{2}$-zero curve passes through the base space.
As illustrated in Fig. \ref{sfig:topo-charge-rel-arg}(b), because
the radius function $\rho\left(\theta_{1}\right)$ must be sandwiched
between $|\beta_{2}^{(\mathrm{a})}|$ and $|\beta_{2}^{(\mathrm{b})}|$,
the direction of the $\beta_{2}$-zero curves is determined by the
derivative of $\left|\beta_{2}^{(\mathrm{a})}\right|$ and $\left|\beta_{2}^{(\mathrm{b})}\right|$
against $\theta_{1}$ at the PMGBZ points. If $\frac{d\left|\beta_{2}^{(\mathrm{a})}\right|^{2}}{d\theta_{1}}|_{p_{\mathrm{a}}}>\frac{d\left|\beta_{2}^{(\mathrm{b})}\right|^{2}}{d\theta_{1}}|_{p_{\mathrm{b}}}$,
$p_{\mathrm{a}}$ takes the positive charge and $p_{\mathrm{b}}$
the negative charge, and vice versa. In general, the derivative of
a $\beta_{2}$-zero is derived from the implicit equation $f\left(E,re^{i\theta_{1}},\beta_{2}^{(u)}\right)=0$,
i.e., 
\begin{equation}
\frac{d\beta_{2}^{(u)}}{d\theta_{1}}=-i\frac{\left.\frac{\partial f}{\partial\beta_{1}}\right|_{p_{u}}re^{i\theta_{1}}}{\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{u}}},\label{seq:derivative-beta-2}
\end{equation}
where $u=\mathrm{a},\mathrm{b}$. Then, the derivative of the absolute
value reads, 
\begin{align}
\frac{d\left|\beta_{2}^{(u)}\right|^{2}}{d\theta_{1}} & =2\text{Re}\left(\beta_{2}^{(u)*}\frac{d\beta_{2}^{(u)}}{d\theta_{1}}\right)=2\text{Im}\left(\frac{\beta_{2}^{(u)*}re^{i\theta_{1}}\left.\frac{\partial f}{\partial\beta_{1}}\right|_{p_{u}}}{\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{u}}}\right).
\end{align}
For simplicity, we denote 
\begin{equation}
K\left(E,\beta_{1},\beta_{2}\right)=\frac{\beta_{2}^{*}\beta_{1}\left.\frac{\partial f}{\partial\beta_{1}}\right|_{\left(E,\beta_{1},\beta_{2}\right)}}{\left.\frac{\partial f}{\partial\beta_{2}}\right|_{\left(E,\beta_{1},\beta_{2}\right)}},\label{eq:K-fun}
\end{equation}
then we have $\frac{d\left|\beta_{2}^{(\mathrm{a})}\right|^{2}}{d\theta_{1}}|_{p_{\mathrm{a}}}=2\text{Im}\left[K\left(E,\mathcal{B}_{1},\mathcal{B}_{2}\right)\right]$
and $\frac{d\left|\beta_{2}^{(\mathrm{b})}\right|^{2}}{d\theta_{1}}|_{p_{\mathrm{b}}}=2\text{Im}\left[K\left(E,\mathcal{B}_{1},\mathcal{B}_{2}e^{i\phi}\right)\right]$.
By comparing the derivatives, the topological charges of $p_{\mathrm{a}}$
and $p_{\mathrm{b}}$ are determined.

Next, we will show the relation between the function $K$ defined
above and the change in the relative phase $\phi$ when the modulus
of $\mathcal{B}_{1}$ increases, i.e., the sign of $\partial\left|\mathcal{B}_{1}\right|^{2}/\partial\phi$
when $E$ is kept constant. The derivative can be calculated from
the auxiliary GBZ equations, 
\begin{align}
f\left(E,\mathcal{B}_{1},\mathcal{B}_{2}\right) & =0,\label{seq:aGBZ-1}\\
f\left(E,\mathcal{B}_{1},\mathcal{B}_{2}e^{i\phi}\right) & =0,\label{seq:aGBZ-2}
\end{align}
as implicit functions of $\mathcal{B}_{1}\left(E,\phi\right)$ and
$\mathcal{B}_{2}\left(E,\phi\right)$. By direct calculations, the
derivative reads, 
\begin{equation}
\frac{\partial\mathcal{B}_{1}}{\partial\phi}=\frac{-ie^{i\phi}\mathcal{B}_{2}\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{a}}}\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{b}}}}{\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{a}}}\left.\frac{\partial f}{\partial\beta_{1}}\right|_{p_{\mathrm{b}}}-e^{i\phi}\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{b}}}\left.\frac{\partial f}{\partial\beta_{1}}\right|_{p_{\mathrm{a}}}},\label{eq:diff-mGBZ-beta1}
\end{equation}
and the derivative of $\left|\mathcal{B}_{1}\right|^{2}$ reads, 
\begin{align}
\frac{\partial\left|\mathcal{B}_{1}\right|^{2}}{\partial\phi} & =2C\text{Im}\left[K\left(E,\mathcal{B}_{1},\mathcal{B}_{2}\right)-K\left(E,\mathcal{B}_{1},\mathcal{B}_{2}e^{i\phi}\right)\right],\nonumber \\
 & =C\left(\left.\frac{d|\beta_{2}^{(\mathrm{a})}|^{2}}{d\theta_{1}}\right|_{p_{\mathrm{a}}}-\left.\frac{d|\beta_{2}^{(\mathrm{b})}|^{2}}{d\theta_{1}}\right|_{p_{\mathrm{b}}}\right),\label{seq:relation-dnorm-K-fun}
\end{align}
where, 
\begin{equation}
C=\frac{\left|\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{b}}}\right|^{2}\cdot\left|\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{a}}}\right|^{2}}{\left|\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{a}}}\left.\frac{\partial f}{\partial\beta_{1}}\right|_{p_{\mathrm{b}}}-e^{i\phi}\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{b}}}\left.\frac{\partial f}{\partial\beta_{1}}\right|_{p_{\mathrm{a}}}\right|^{2}},
\end{equation}
is a positive value. Therefore, the sign of $\partial\left|\mathcal{B}_{1}\right|^{2}/\partial\phi$
is the same as the topological charge of $p_{\mathrm{a}}$, and opposite
to the topological charge of $p_{\mathrm{b}}$.

For the strip winding number $W_{\text{strip}}\left(E,r\right)$,
the sign $\left(-1\right)^{\tau_{j}}$ equals $+1$ when $\phi_{j}$
starts from a positive charge, and equals $-1$ when starting from
a negative charge. According to the discussion above, $\left(-1\right)^{\tau_{j}}$
equals the sign of $\partial\left|\mathcal{B}_{1}\right|^{2}/\partial\phi$.
That is, when the radius $r$ increases to some $r^{\prime}>r$, the
sign of $\phi_{j}^{\prime}-\phi_{j}$ equals $\left(-1\right)^{\tau_{j}}$.
As a result, Eq. (\ref{seq:W-strip}) is simplified into, 
\begin{equation}
W_{\text{strip}}\left(E,r^{\prime}\right)-W_{\text{strip}}\left(E,r\right)=\frac{1}{2\pi}\sum_{j=1}^{k}\left|\phi_{j}^{\prime}-\phi_{j}\right|,
\end{equation}
which proves the claim in the main text.

\section{Relation of the Degrees in Quasi-1D and 2D Characteristic Polynomials}

\label{app:relation-degrees}

In the main text, we state without proof that the lowest and highest
degrees of $\beta_{1}$ in the quasi-1D ChP $F(E,\beta_{1})$ and
2D ChP $f(E,\beta_{1},\beta_{2})$ are related to each other by, 
\begin{align}
M & =M_{1}L_{2}+O(1),\label{seq:relation-M}\\
N & =N_{1}L_{2}+O(1).\label{seq:relation-N}
\end{align}
In this section, we will prove the two relations for a general 2D
lattice.

First, we expand the 2D non-Bloch Hamiltonian $h(\beta_{1},\beta_{2})$
as the polynomial of $\beta_{2}$, i.e., 
\begin{equation}
h\left(\beta_{1},\beta_{2}\right)=\sum_{j=-M_{2}}^{N_{2}}h^{(j)}\left(\beta_{1}\right)\beta_{2}^{j},\label{seq:h-in-polynomial}
\end{equation}
where $h^{(j)}(\beta_{1})$ is an $m\times m$ matrix. Then, the hybrid
non-Bloch Hamiltonian can be expanded into the block Toeplitz form
under the single particle basis, i.e., 
\begin{equation}
\mathcal{H}\left(\beta_{1}\right)=\underbrace{\left(\begin{matrix}h^{(0)} & h^{(-1)} & \cdots & h^{(-M_{2})}\\
h^{(1)} & h^{(0)} & h^{(-1)} & \cdots & \ddots\\
\vdots & h^{(1)} & h^{(0)} & h^{(-1)} & \cdots & h^{(-M_{2})}\\
h^{(N_{2})} & \vdots & h^{(1)} & h^{(0)} & \ddots & \vdots\\
 & \ddots & \vdots & \ddots & \ddots & h^{(-1)}\\
 &  & h^{(N_{2})} & \cdots & h^{(1)} & h^{(0)}
\end{matrix}\right)}_{L_{2}\text{ blocks}}.\label{seq:quasi-1D-Toeplitz}
\end{equation}

We first consider the lowest degrees. The case of the highest degrees
can be proved by the same reasoning. It is noted that the matrix elements
of $h^{(j)}(\beta_{1})$ are all Laurent polynomials of $\beta_{1}$.
Assume that the lowest degree of $\beta_{1}$ in $h_{\mu\nu}(\beta_{1},\beta_{2})$
is $-M_{\mu\nu}$, and the degree of $\beta_{2}$ in the term with
lowest degree of $\beta_{1}$ is $t_{\mu\nu}$, then, the elements
of the ChP $f(E,\beta_{1},\beta_{2})$ can be expressed as the leading
term $a_{ij}\beta_{1}^{-M_{ij}}\beta_{2}^{t_{ij}}$ plus a remainder
in which the degrees of $\beta_{1}$ are greater than $-M_{ij}$,
like Eq. (\ref{seq:expansion-f}), 
\begin{align}
f\left(E,\beta_{1},\beta_{2}\right) & =\det\left[E-h\left(\beta_{1},\beta_{2}\right)\right],\nonumber \\
 & =\begin{vmatrix}a_{11}\beta_{1}^{-M_{11}}\beta_{2}^{t_{11}}+\cdots & a_{12}\beta_{1}^{-M_{12}}\beta_{2}^{t_{12}}+\cdots & \cdots & a_{1m}\beta_{1}^{-M_{1m}}\beta_{2}^{t_{1m}}+\cdots\\
a_{21}\beta_{1}^{-M_{21}}\beta_{2}^{t_{21}}+\cdots & a_{22}\beta_{1}^{-M_{22}}\beta_{2}^{t_{22}}+\cdots &  & \vdots\\
\vdots &  & \ddots\\
a_{m1}\beta_{1}^{-M_{m1}}\beta_{2}^{t_{m1}}+\cdots & \cdots &  & a_{mm}\beta_{1}^{-M_{mm}}\beta_{2}^{t_{mm}}+\cdots
\end{vmatrix},\label{seq:expansion-f}
\end{align}
where $a_{ij}\in\mathbb{C}$ are the coefficients. In the determinant
of Eq. (\ref{seq:expansion-f}), only the leading terms contribute
to the lowest degree of $f(E,\beta_{1},\beta_{2})$. Therefore, the
term in $f(E,\beta_{1},\beta_{2})$ with lowest degree in $\beta_{1}$
can be construct by the following steps: 
\begin{enumerate}
\item Find $i_{\max},j_{\max}$ that maximizes $M_{ij}$, i.e. $M_{i_{\max}j_{\max}}=\max_{i,j}\left\{ M_{ij}\right\} $. 
\item Get the $(i_{\max},j_{\max})$-cofactor of $f\left(E,\beta_{1},\beta_{2}\right)$,
denoted as $f_{1}(E,\beta_{1},\beta_{2})$, then factorize $f(E,\beta_{1},\beta_{2})$
as, 
\begin{equation}
f\left(E,\beta_{1},\beta_{2}\right)=a_{i_{\max}j_{\max}}\beta_{1}^{-M_{i_{\max}j_{\max}}}\beta_{2}^{t_{i_{\max}j_{\max}}}f_{1}\left(E,\beta_{1},\beta_{2}\right)+\dots
\end{equation}
\item Substitute the determinant $f$ by the cofactor $f_{1}$, and repeat
the two steps above to get $f_{2}$, $f_{3}$, $\dots$, until the
order of the cofactor is reduced to $1$. 
\end{enumerate}
After this procedure, we get a series of indices $(i_{\max}^{(1)},j_{\max}^{(1)}),(i_{\max}^{(2)},j_{\max}^{(2)}),\dots,(i_{\max}^{(m)},j_{\max}^{(m)})$.
According to our construction, $i_{\max}^{(1)},i_{\max}^{(2)},\dots,i_{\max}^{(m)}$
and $j_{\max}^{(1)},j_{\max}^{(2)},\dots,j_{\max}^{(m)}$ are permutations
of $1,2,\dots,m$. Therefore, we can sort $M_{i_{\max}^{(1)}j_{\max}^{(1)}},M_{i_{\max}^{(2)}j_{\max}^{(2)}},\dots,M_{i_{\max}^{(m)}j_{\max}^{(m)}}$
to $M_{1\nu_{1}},M_{2\nu_{2}},\dots,M_{m\nu_{m}}$, and ensure that
$\nu_{1},\nu_{2},\dots,\nu_{m}$ is a permutation of $1,2,\dots,m$.
Then, the minimum negative degree of $f(E,\beta_{1},\beta_{2})$ reads,
\begin{equation}
M_{1}=\sum_{j=1}^{m}M_{j\nu_{j}}.
\end{equation}

Back to the hybrid non-Bloch Hamiltonian $\mathcal{H}(\beta_{1})$
of Eq. (\ref{seq:quasi-1D-Toeplitz}), according to Eq. (\ref{seq:h-in-polynomial}),
the term $a_{j\nu_{j}}\beta_{1}^{-M_{j\nu_{j}}}$ is available in
the matrix element $h_{j\nu_{j}}^{(t_{j\nu_{j}})}(\beta_{1})$. Therefore,
in each row of the blocks in Eq. (\ref{seq:quasi-1D-Toeplitz}), we
can always pick the terms $a_{1\nu_{1}}\beta_{1}^{-M_{1\nu_{1}}},a_{2\nu_{2}}\beta_{2}^{-M_{2\nu_{2}}},\dots,a_{m\nu_{m}}\beta_{2}^{-M_{m\nu_{m}}}$,
except for the first $N_{2}$ or last $M_{2}$ rows. For different
rows of the blocks, the column indices of the selected terms do not
coincide with each other, in that $\nu_{1},\nu_{2},\dots,\nu_{m}$
is a permutation of $1,2,\dots,m$. Next, considering the determinant
$F(E,\beta_{1})\equiv\det[E-\mathcal{H}(\beta_{1})]$, except for
the first $N_{2}$ and last $M_{2}$ rows of the blocks, the lowest
degree of $\beta_{1}$ contributed from each row reads, 
\[
\beta_{1}^{-M_{1\nu_{1}}}\beta_{1}^{-M_{2\nu_{2}}}\cdots\beta_{1}^{-M_{m\nu_{m}}}=\beta_{1}^{-M_{1}},
\]
and the total contribution of the $L_{2}-M_{2}-N_{2}$ rows is, 
\[
\beta_{1}^{-M_{1}\left(L_{2}-M_{2}-N_{2}\right)}=\beta_{1}^{-M_{1}L_{2}+O(1)}.
\]
Because the terms of $\beta_{1}$ in the first $N_{2}$ and last $M_{2}$
rows can only change the total degree of $\beta_{1}$ by a finite
amount (i.e., independent of $L_{2}$), the lowest degree of $F(E,\beta_{1})$
satisfies Eq. (\ref{seq:relation-M}). By the same reasoning, we can
also prove Eq. (\ref{seq:relation-N}).

\section{Numerical Calculation of Quasi-1D Winding Number and Strip Winding
Number}

\label{app:num-W-strip}

\begin{figure}
\centering

\includegraphics{Figures/winding-number-numerical-20250212}

\caption{Comparison of $W(E,r)$ and $W_{\text{strip}}(E,r)$ in 2D HN model
and non-Hermitian Haldane model. (a) Sketch of 2D HN model, where
the parameters are $J_{x1}=1+i$, $J_{x2}=1.5+1.2i$, $J_{y1}=-1+i$
and $J_{y2}=-1.2-0.5i$. The major (minor) axis is marked by the lattice
vector $\mathbf{a}_{1}$ ($\mathbf{a}_{2}$). (b) Numerical calculation
of $W(E,r)/L_{2}$ (blue solid lines) and $W_{\text{strip}}(E,r)$
(orange dashed lines) of 2D HN model with major axis $\mathbf{a}_{1}$.
The width is $L_{2}=40$, and the reference energies $E$ are (i)
$1.00296-0.21641i$, (ii) $1.55832+0.91741i$, (iii) $-2.57608+1.13451i$
and (iv) $-1.57168+0.02125i$. (c) Sketch of non-Hermitian Haldane
model, where the parameters are $t_{1}=0.70502$, $t_{2}=-1.32760$,
$\gamma=2.15618$, $\varphi=0.05877$ and $m=-0.64569$. The unit
cell is marked by the purple ellipse, and the major (minor) axis is
marked by the lattice vector $\mathbf{a}_{1}$ ($\mathbf{a}_{2}$).
(d) Numerical calculation of $W(E,r)/L_{2}$ (blue solid lines) and
$W_{\text{strip}}(E,r)$ (orange dashed lines) of the non-Hermitian
Haldane model with major axis $\mathbf{a}_{1}$. The width is $L_{2}=20$
($40$ sites in the supercell), and the reference energies $E$ are
(i) $0.01162+0.68736i$, (ii) $0.54100-1.99811i$, (iii) $0.59046+1.89903i$
and (iv) $0.10591-0.34594i$.}
\label{sfig:numerical-comparison-winding} 
\end{figure}

In the main text, we have derived the relation between the quasi-1D
winding number $W(E,r)$ and the strip winding number $W_{\text{strip}}(E,r)$,
which reads $W_{\text{strip}}(E,r)=W(E,r)/L_{2}+O(1/L_{2})$. In order
to illustrate the relation, we numerically calculate $W(E,r)$ and
$W_{\text{strip}}(E,r)$ in two specific models: the 2D HN model and
the non-Hermitian Haldane model, and compare the curves of $W(E,r)/L_{2}$
and $W_{\text{strip}}(E,r)$ for randomly selected reference energies.

Figure \ref{sfig:numerical-comparison-winding}(a) shows the 2D HN
model with complex coupling coefficients, the same as the example
in the main text. In the numerical calculation, the coefficients are
the same as the numerical examples of the main text, which are $J_{x1}=1+i$,
$J_{x2}=1.5+1.2i$, $J_{y1}=-1+i$ and $J_{y2}=-1.2-0.5i$. The lattice
vector $\mathbf{a}_{1}=(1,1)$ is selected as the main axis and $\mathbf{a}_{2}=(0,1)$
as the minor axis. To calculate the quasi-1D winding number, the supercell
is selected by selecting successive $L_{2}$ unit cells along the
minor axis $\mathbf{a}_{2}$. Figure \ref{sfig:numerical-comparison-winding}(b)
illustrates the numerical result of $W_{\text{strip}}(E,r)$ and $W(E,r)/L_{2}$
of the 2D HN model with randomly generated reference energies. In
the numerical calculation, $L_{2}$ is set to $40$, and the reference
energies are (i) $1.00296-0.21641i$, (ii) $1.55832+0.91741i$, (iii)
$-2.57608+1.13451i$ and (iv) $-1.57168+0.02125i$. As shown in Fig.
\ref{sfig:numerical-comparison-winding}(b), the curves of $W_{\text{strip}}(E,r)$
against $r$ (orange dashed curves) are piecewise smooth curves, while
the curves of $W(E,r)/L_{2}$ (blue solid curves) show platforms because
$L_{2}$ is finite. For all four samples in Fig. \ref{sfig:numerical-comparison-winding}(b),
$W(E,r)/L_{2}$ fits well with $W_{\text{strip}}(E,r)$.

To show the generality of our conclusion, we also compare $W(E,r)$
and $W_{\text{strip}}(E,r)$ in a non-Hermitian version of the Haldane
model. As shown in Fig. \ref{sfig:numerical-comparison-winding}(c),
the nearest coupling is $t_{1}\in\mathbb{R}$, and the next-nearest
coupling has a nonreciprocal phase, i.e., the coupling coefficient
is $t_{2}e^{i\left(\varphi+\gamma\right)}$ along the direction of
the arrow, and $t_{2}e^{-i\left(\varphi-\gamma\right)}$ against the
arrow. Under the basis $\{\mathbf{a}_{1},\mathbf{a}_{2}\}$ shown
in Fig. \ref{sfig:numerical-comparison-winding}(c), the Hamiltonian
reads, 
\begin{align}
H= & m\sum_{r_{1},r_{2}}\left(c_{r_{1},r_{2},\mathrm{a}}^{\dagger}c_{r_{1},r_{2},\mathrm{a}}-c_{r_{1},r_{2},\mathrm{b}}^{\dagger}c_{r_{1},r_{2},\mathrm{b}}\right)+\nonumber \\
 & +t_{1}\sum_{r_{1},r_{2}}\left(c_{r_{1},r_{2},\mathrm{b}}^{\dagger}c_{r_{1},r_{2},\mathrm{a}}+c_{r_{1},r_{2}+1,\mathrm{b}}^{\dagger}c_{r_{1},r_{2},\mathrm{a}}+c_{r_{1}-1,r_{2},\mathrm{b}}^{\dagger}c_{r_{1},r_{2},\mathrm{a}}+\text{h.c.}\right)+\nonumber \\
 & +t_{2}e^{i\gamma}\sum_{r_{1},r_{2}}\left[e^{i\varphi}\left(c_{r_{1}-1,r_{2},\mathrm{a}}^{\dagger}c_{r_{1},r_{2},\mathrm{a}}+c_{r_{1},r_{2}-1,\mathrm{a}}^{\dagger}c_{r_{1},r_{2},\mathrm{a}}+c_{r_{1}+1,r_{2}+1,\mathrm{a}}^{\dagger}c_{r_{1},r_{2},\mathrm{a}}\right)+\right.\nonumber \\
 & \left.+e^{-i\varphi}\left(c_{r_{1}-1,r_{2},\mathrm{b}}^{\dagger}c_{r_{1},r_{2},\mathrm{b}}+c_{r_{1},r_{2}-1,\mathrm{b}}^{\dagger}c_{r_{1},r_{2},\mathrm{b}}+c_{r_{1}+1,r_{2}+1,\mathrm{b}}^{\dagger}c_{r_{1},r_{2},\mathrm{b}}\right)+\text{h.c.}\right],
\end{align}
where $c_{r_{1},r_{2},\mu},r_{1}\in\mathbb{Z},r_{2}\in\mathbb{Z},\mu=\mathrm{a},\mathrm{b}$
is the annihilator at the sublattice $\mu$ in the unit cell with
coordinate $r_{1}\mathbf{a}_{1}+r_{2}\mathbf{a}_{2}$, and $m\in\mathbb{R}$
is the detuning between the site $\mathrm{a}$ and the site $\mathrm{b}$.
In the numerical calculation, the parameters $t_{1},t_{2},\gamma,\varphi$
and $m$ are randomly generated as $t_{1}=0.70502$, $t_{2}=-1.32760$,
$\gamma=2.15618$, $\varphi=0.05877$, and $m=-0.64569$. As shown
in Fig. \ref{sfig:numerical-comparison-winding}(d), we calculate
$W(E,r)$ and $W_{\text{strip}}(E,r)$ against $r$ in the direction
of $\mathbf{a}_{1}$ for four random reference energies, which are
(i) $0.01162+0.68736i$, (ii) $0.54100-1.99811i$, (iii) $0.59046+1.89903i$
and (iv) $0.10591-0.34594i$. For the quasi-1D winding number $W(E,r)$,
$L_{2}$ is set to $20$ (that is, $40$ sites in the supercell).
Similarly to the case of the 2D HN model, the curves of $W(E,r)/L_{2}$
fit well with $W_{\text{strip}}(E,r)$ for the four samples.

\section{Amoeba spectrum of the 2D HN model}

\label{app:calculation-Amoeba}

In this section, we will prove that the Amoeba spectrum of the 2D
HN model is the same as the $x$($y$)-SGBZ spectrum, i.e., 
\begin{equation}
\sigma_{\text{Amoeba}}=\left\{ E\in\mathbb{C}\mid E=2e^{i\delta_{x}}\text{Re}\left(J_{x}^{*}e^{i\theta_{x}}\right)+2e^{i\delta_{y}}\text{Re}\left(J_{y}^{*}e^{i\theta_{y}}\right),\theta_{x},\theta_{y}\in\left[-\pi,\pi\right]\right\} ,\label{seq:Amoeba-spectrum}
\end{equation}
shown as Fig. \ref{sfig:calculation-Amoeba}(a). To prove Eq. (\ref{seq:Amoeba-spectrum}),
we need to prove that Amoeba $\mathcal{A}(E)$ exhibits a central
hole when $E$ is outside the spectrum and has no central holes when
$E$ is inside the spectrum.

When $|\beta_{x}|=\exp(\mu_{x})$ and $|\beta_{y}|=\exp(\mu_{y})$,
the ChP reads, 
\begin{align}
f_{xy}\left(E,e^{\mu_{x}+i\theta_{x}},e^{\mu_{y}+i\theta_{y}}\right) & =E-e^{i\delta_{x}}\left(\gamma_{x}e^{-\mu_{x}}J_{x}e^{-i\theta_{x}}+\gamma_{x}^{-1}e^{\mu_{x}}J_{x}^{*}e^{i\theta_{x}}\right)-e^{i\delta_{y}}\left(\gamma_{y}e^{-\mu_{y}}J_{y}e^{-i\theta_{y}}+\gamma_{y}^{-1}e^{\mu_{y}}J_{y}^{*}e^{i\theta_{y}}\right),\nonumber \\
 & =E-U_{x}\left(\theta_{x}\right)-U_{y}\left(\theta_{y}\right),
\end{align}
where the functions $U_{j}(\theta_{j}),j=x,y$ are defined as, 
\begin{align}
U_{j}\left(\theta_{j}\right) & \equiv e^{i\delta_{j}}\left(\gamma_{j}e^{-\mu_{j}}J_{j}e^{-i\theta_{j}}+\gamma_{j}^{-1}e^{\mu_{j}}J_{j}^{*}e^{i\theta_{j}}\right),\nonumber \\
 & =e^{i\delta_{j}}\left[\left(\gamma_{j}e^{-\mu_{j}}+\gamma_{j}^{-1}e^{\mu_{j}}\right)\text{Re}\left(J_{j}^{*}e^{i\theta_{j}}\right)+i\left(\gamma_{j}^{-1}e^{\mu_{j}}-\gamma_{j}e^{-\mu_{j}}\right)\text{Im}\left(J_{j}^{*}e^{i\theta_{j}}\right)\right].
\end{align}
The trajectory of $v_{j}(\theta_{j})$ is an ellipse on the complex
plane, whose major axis is parallel to $\exp(i\delta_{j})$.

By definition, $(\mu_{x},\mu_{y})$ belongs to $\mathcal{A}(E)$ when
and only when the ChP $f_{xy}(E,e^{\mu_{x}+i\theta_{x}},e^{\mu_{y}+i\theta_{y}})$
vanishes for some $\theta_{x}$ and $\theta_{y}$, or equivalently
when the ellipse $E-U_{x}(\theta_{x})$ intersects with the ellipse
$U_{y}(\theta_{y})$. Figure \ref{sfig:calculation-Amoeba}(b) and
\ref{sfig:calculation-Amoeba}(c) shows the cases for $E_{1}\notin\sigma_{\text{Amoeba}}$
and $E_{2}\in\sigma_{\text{Amoeba}}$, respectively. According to
the geometric-arithmetic mean inequality, the length of the major
semi-axis is $(\gamma_{j}e^{-\mu_{j}}+\gamma_{j}^{-1}e^{\mu_{j}})|J_{j}|\geq2|J_{j}|$.
Therefore, the ellipse $U_{j}(\theta_{j})$ must enclose (or coincide
with) the segment $2e^{i\delta_{j}}\text{Re}(J_{j}^{*}e^{i\theta_{j}}),\theta_{j}\in[-\pi,\pi]$,
shown as the dotted segments in Fig. \ref{sfig:calculation-Amoeba}(b)
and \ref{sfig:calculation-Amoeba}(c).

For $E_{1}\notin\sigma_{\text{Amoeba}}$, as shown in Fig. \ref{sfig:calculation-Amoeba}(b),
the two dotted segments have no intersections. When $(\mu_{x},\mu_{y})$
is close to $(\ln\gamma_{x},\ln\gamma_{y})$ enough, the ellipses
of $E_{1}-U_{x}(\theta_{x})$ and $U_{y}(\theta_{y})$ are separated
to each other, so that $(\mu_{x},\mu_{y})\notin\mathcal{A}(E_{1})$
in this case. For the derivative of the Ronkin functions, both the
winding number of $E_{1}-U_{x}(\theta_{x})$ around $U_{y}(\theta_{y})$,
and the winding number of $U_{y}(\theta_{y})$ around $E_{1}-U_{x}(\theta_{x})$,
vanish for arbitrary $\theta_{x}$ and $\theta_{y}$, so that the
derivative of the Ronkin function vanishes in this case. Therefore,
the Amoeba $\mathcal{A}(E_{1})$ shows a central hole in the neighborhood
of $(\ln\gamma_{x},\ln\gamma_{y})$.

For $E_{2}\in\sigma_{\text{Amoeba}}$, as shown in Fig. \ref{sfig:calculation-Amoeba}(c),
the two dotted segments intersect each other. Therefore, the two ellipses
of $E_{2}-U_{x}(\theta_{x})$ and $U_{y}(\theta_{y})$ do not have
intersections only when one ellipse encloses the other. When the two
ellipses intersect, $(\mu_{x},\mu_{y})\in\mathcal{A}(E_{2})$. When
one ellipse encloses the other, the winding number of the outer ellipse
around the points on the inner ellipse is non-zero, resulting in the
non-zero gradient of the Ronkin function. Consequently, $\mathcal{A}(E_{2})$
does not have central holes.

\begin{figure}
\centering

\includegraphics{Figures/sfig-calculation-Amoeba}

\caption{Sketches for the Amoeba spectrum of the 2D HN model. (a) Amoeba spectrum
of the 2D HN model. $E_{1}$ and $E_{2}$ are two reference energies
out of and in the spectrum, respectively. (b, c) Sketches of the ellipses
$E-U_{x}(\theta_{x})$ and $U_{y}(\theta_{y})$ for (b) $E=E_{1}$
and (c) $E=E_{2}$. The black dotted segments enclosed by the green
and orange ellipses denote the trajectories of $E-2e^{i\delta_{x}}\text{Re}(J_{x}^{*}e^{i\theta_{x}})$
and $2e^{i\delta_{y}}\text{Re}(J_{y}^{*}e^{i\theta_{y}})$, respectively.}
\label{sfig:calculation-Amoeba} 
\end{figure}

\end{document}

